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Preface 


The  first  step  in.  problem  solving  is  determining  the  reasons  and 
causes  that  underlie  the  problem.  Such,  has  been  the  purpose  of  this  dis¬ 
sertation.  Previous  work  involving  the  second  order  even-parity  form  of 
the  Boltzmann  transport  equation  successfully  demonstrated  its  ability 
to  generate  ray  effect  free  solutions.  These  solutions  were  obtained  in 
transport  media  requiring  relatively  simple  geometry  and  isotropic 
scatter.  The  obvious  extension  of  this  work  was  to  attempt  a  more  complex 

*  o 

transport  problem.  An  air-over-ground  problem  was  chosen  since  it  adds 
the  necessary  complication  to  meaningfully  test  the  feasibility  of  such 
an  extension.  Additionally,  the  solution  of  this  problem  is  important 
in  predicting  the  survivability  and  vulnerability  of  military  systems 
in  a  nuclear  environment.  Thus,  a  successful  extension  would  yield  a 
worthwhile  product.  Unfortunately  such  an  extension  was  not  possible 
using  the  standard  solution  techniques  currently  available.  The  positive 
aspect  of  this  research  is  that  the  problems  associated  with  applying 
this  equation  to  more  complex  transport  problems  are  identified  and  the 
causes  are  known.  The  results  of  this  research  will  provide  direction 
for  future  efforts  aimed  at  developing  new  solution  techniques  that  will 
successfully  solve  the  even-parity  form  of  the  Boltzmann  equation  formu¬ 
lated  for  these  more  complex  problems. 

As  in  all  research  efforts,  several  people  played  an  important  role. 
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cant.  I  would  also  like  to  express  my  gratitude  to  Dr.  Bridgman,  Dr.  Jones 
and  Dr.  Kelleher,  the  other  members  of  my  committee,  for  their  suggestions 
and  constant  encouragement.  Two  members  of  the  Air  Force  Weapons  Labora¬ 
tories  (AFWL)  deserve  special  recognition,  Mr.  John  Burgio  and  Mr.  Harry 
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Abstract 

Published  work  indicated  that  a  finite  element  solution  of  the  even- 
parity  form  of  the  Boltzmann  equation  (.EPFBE)  provided  a  means  for  curing 
the  ray  effect  in  transport  problems.  This  conclusion  was  made  after 
examining  the  solutions  of  some  simply  defined  problems  involving  plane 
geometry  and  isotropic  scatter.  The  purpose  of  this  research  was  to 
determine  the  feasibility  of  applying  this  equation  to  more  complex  prob¬ 
lems.  The  air-over-ground  problem  was  chosen.  This  problem  is  important 

*  * 

in  nuclear  weapon  effects  calculations  and  requires  two-dimensional 
cylindrical  geometry,  anisotropic  scatter,  and  a  multigroup  solution. 
Additionally,  the  discrete  ordinates  method  applied  to  this  problem  gener¬ 
ates  solutions  with  severe  ray  effects.  Several  numerical  approaches 
based  on  the  finite  element  method  were  attempted.  The  Galerkin  method 
was  applied  to  the  weak  form  of  the  EPFBE.  A  bilinear  Lagrange  polynomial 
trial  solution  and  specially  derived  synthesis  function  trial  solution 
were  used.  The  Galerkin  method  proved  infeasible  because  the  integrals 
resulting  from  this  method  could  not  be  efficiently  evaluated  either 
numerically  or  analytically.  To  solve  this  integration  problem,  the 
collocation  method  was  attempted.  A  trial  solution  consisting  of  cubic 
splines  and  a  simply  defined  angular  synthesis  function  was  used.  The 
collocation  method  allowed  the  analytic  evaluation  of  the  resulting  inte¬ 
grals  but  forced  a  fixed  anisotropy  on  the  solution.  The  multi  group 
method  applied  to  the  EPFBE  resulted  in  a  nested  integral  problem  involv¬ 
ing  the  source  terms  of  this  equation.  The  complexity  of  this  nesting 
problem  increased  proportionately  with  the  number  of  energy  group  used. 
This  research  demonstrated  that  the  finite  element  method  cannot  be  cost 
effectively  used  in  solving  the  EPFBE  for  transport  problems  requiring 
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complex  geometries,  anisotropic  scatter,  and  a  multigroup  solution. 
Criteria  were  developed  from  this  research  that  provides  guidelines  for 
pursuing  future  work  related  to  the  EPFBE.  Recommendations  based  on 
these  criteria  are  made. 
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T.  Introduction 


The.  steady-state,  transport  of  neutrons  and  gamma  rays  from  a  point 
source  of  simulated  nuclear  weapon  radiat.on  in  the  atmosphere  is  a  prob¬ 
lem  of  fundamental  importance  in  nuclear  effects  calculations.  The  equa¬ 
tion  whose  solution  prescribes  the  steady-state,  neutral -particle  distri¬ 
bution  in  the  phase  space  of  position,  particle  direction,  and  particle 
energy  is  the  Boltzmann  transport  equation.  This  equation  is  a  linear 

form  of  the  Boltzmann  equation  used  in  several  fields  of  physics  and  is 

* 

a  statement  of  particle  conservation  as  applied  to  an  infinitesimal  ele¬ 
ment  of  volume,  direction  and  energy,  and  may  be  written  as 


-a.-  v  "Kcf.A.e)  +  crT(r,ej  3[c?,sL,e> 
where 


=  ScrA.,E)  +• 

^ d  A C£ (r,  E)  r,  A,'£ ') 


(l.D 


r  =  the  spatial  position  vector,  in  units  of  distance, 

£  =  the  particle  energy, 

A.  =  a  unit  vector  in  the  direction  of  particle  motion, 

V  =  the  gradient  operator 

S(r,A.,el  ~  The  particle  source  density,  in  particles/unit  volume/unit 
energy, 

^cr.H.e)  ~  The  angular  particle  fluence,  in  particles/unit  area/ 
steradian/unit  energy, 

crT(r,El  =  the  macroscopic  total  interaction  cross  section  at  posi¬ 
tion  r  and  for  energy  E,  per  unit  distance,  per  unit 
energy, 

=  the  macroscopic  differential  scattering  cross  section 
at  position  r  for  a  particle  of  direction  A'  with  energy 
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/  A 

E  scattering  to  a  di  recti  on  J\.  and  energy  E.  The  units 
are  per  unit  distance  per  steradian  per  unit  energy. 

Analytical  solutions  to  this  equation  occur  in  only  special  cases 
where  simplifying  assumptions  can  be  applied  (Reference  1).  Unfortunately, 
these  assumptions  do  not  generate  a  realistic  scenario  for  performing 
nuclear  effects  calculations,  thus  requiring  the  use  of  numerical  methods 
in  these  types  of  transport  problems. 

The  scenario  most  used  in  performing  nuclear  effects  calculations 

*  * 

is  a  nuclear  air  burst  detonation.  Analytically  this  scenario  is  modeled 
by  a  two-dimensional  cylindrical  geometry  (azimuthal  symmetry  assumed) 
composed  of  air  with  ah  exponentially  varying  density  over  a  flat  ground. 
This  type  of  configuration  is  referred  to  as  air-over-ground  geometry. 

In  this  type  of  configuration,  there  is  an  interface  at  the  air-ground 
boundary  across  which  a  discontinuous  change  in  density  occurs.  The 
importance  of  including  the  ground  was  demonstrated  by  Straker  (Ref.  2) 
who  showed  a  significant  effect  on  the  atmospheric  neutron  distribution 
due  to  the  presence  of  this  relatively  high  density  medium. 

The  conservative  form  of  the  Boltzmann  transport  equation  for  the 
selected  cylindrical  geometry  can  be  written  as 

J  coscx)  d(e Hce.i.Z.E))  _  _L 

e  df  e  J* 

+  m  $  •+•  <rT(e,i,E)  ice.i.A.Ei  ~ 


'Sie,*,A,e)  +)  Je' I  a'-»A .£*-»£) P,*, A\e') 

■i  An 


(1.2) 


where 

6  =  the  distance  in  the  radial  spatial  direction 
2  =  the  distance  in  the  axial  spatial  direction 
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JLK  =  direction  cosine  with  respect  to  the  Z-axis 
OC  =  azimuthal  angle  describing  particle  direction 
Figure  1  describes  the  coordinate  relation  between p ,  1  ,*t,X . 

In  Eq.  (1.2)  the  energy  dependence  of  the  angular  particle  fluence 
is  represented  as  a  continuous  function.  In  practice  the  energy  variable 
is  divided  into  a  number  of  finite*  discrete  energy  groups.  This  discre¬ 
tization  of  the  energy  variable  is  known  as  the  multigroup  method  (Ref.  1) 
and  is  used  to  derive  the  set  of  multigroup  transport  equations  by  inter¬ 
grating  Eq.  (1.2)  over  each  discrete  energy  group.  Applying  the  muTtigroup 
method  to  Eq.  (1.2)  yields 

>1  7 -  u* '  co s(x)  d(e  «.  JL  4  ul  ’ Tfry.n)  L 

e  e  dx 


+ 


u. 


6* 


5-1 


(1.3) 


where 

9  =  the  energy  group  designator,  g=l  being  highest  energy  group 

G  =  total  number  of  energy  groups 

=  the  angular  particle  fluence  for  energy  group  g,  and 
*  the  macroscopic  differential  scattering  cross  section 
at  position  (p,  z)  for  a  particle  of  direction  (/t*,  x) 
in  energy  group  g’  scattering  to  a  direction  ((4,  X)  and 
to  energy  group  g. 

To  complete  the  specification  of  this  analytical  problem,  boundary 
conditions  must  be  assigned  to  Eq.  (1.3).  One  boundary  condition  often 
used  is  referred  to  as  a  vacuum  boundary  condition  and  is  represented  as 
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Fig  1.  Coordinate  Systems  Used  In  Neutral  Particle  Transport 
Problems  For  2-D  Cylindrical  Geometry 


J  :  O  ,  A-  a.  <  O 


(1.4) 


where 

-A. 

Ksv  =  the  position  vector  representing  a  spatial  location  on  an 
exterior  surface 

ri  =  a  unit  vector  normal  to  the  exterior  surface  located  at 
The  second  boundary  condition  is  applied  along  the  axis  of  the  cylinder 
at  £3=0  and  is  referred  to  as  the  symmetry  or  reflective  boundary  condi¬ 
tion.  This  boundary  condition  can  be  written  as 


3T  A)  =  I  (O, 

where  is  defined  by  a-  a'  =-a-A  ,  and(A*>r;- a  =  O  . 


(1.5) 


A  nuclear  burst  is  represented  as  an  isotropic  point  source  in  the 
atmosphere.  In  the  cylindrical  air-over-ground  geometry  the  isotropic 
point  source  is  located  along  the  axis  of  the  cylinder  (p=0)  at  the 
desired  burst  altitude  (Zg).  Mathematically  this  point  source  can  be 
represented  as 

-So*  cTcfJ 

H  7i  (1*6) 

where 

«S0  =  isotropic  source  in  group  g 
S(i-2,)Y$(e)  =  Dirac  delta  function 

During  the  course  of  this  dissertation  it  became  necessary  to  define 
another  source  term  known  as  the  first-scatter  source.  This  source  can 


be  expressed  as 


I  ■»  -*■  I  * 


Hr, 


(1.7) 


where 
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rs  -  position  vector  of  the  nuclear  burst 
t"  -  position  vector  In  spatial  domain 

A.e  =  direction  vector  aligned  with  ther-£  vector  at  the  spatial 
position  r\ 

A  =  /  O'* ds  along  the  ray  defined  byAe» 

See  Fig.  2  for  a  clarification  of  the  above  definitions.  The  reason  for 

using  the  first-scatter  source  will  be  explained  in  latter  chapters. 

Combining  the  previously  derived  results  leads  to  the  following 

«  • 

form  of  the  Boltzmann  transport  equation: 


The  sum  over  energy  groups  has  been  terminated  at  g  because  scatter 
occurs  only  from  higher  to  equal  or  lower  energies. 

The  most  common  techniques  used  to  solve  the  Boltzmann  equation  for 
the  air-over-ground  problem  are  Monte  Carlo,  discrete  ordinates,  and  Mass 
Integral  Scaling  (MIS).  Each  method  is  deficient  in  some  respect.  Monte 
Carlo  demands  a  substantial  amount  of  computational  time  due  to  its 
statistical  nature  (Ref.  3).  Discrete  ordinates  (Refs.  4  and  5)  is 
susceptible  to  a  computational  anomaly  called  the  ray  effect.  All  attempts 
to  eliminate  or  mitigate  this  computation  anomaly  in  the  discrete  ordi¬ 
nates  method  have  resulted  in  a  substantial  increase  in  computational 
time  (Refs.  6-12).  MIS  (Refs.  13-15)  cannot  provide  accurate  results 


&D>  STKEAMMG  MtECTVOM  AT  (f,Z) 
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Fig  2.  First  Scatter  Source  Coordinate  System  Definition 


near  the  air-ground  interface  at  slant  ranges  less  than  1000  feet.  The 
deficiencies  exhibited  by  these  commonly  used  solution  techniques, 
motivated  several  new  approaches  for  solving  transport  problems. 
Even-Parity  Equation  and  Space  Angle  Synthesis 

In  1961  Vladimirov  (Ref.  16)  derived  a  self-adjoint  form  of  the 
Boltzmann  transport  equation  known  as  the  second-order  even-parity  flu- 
ence  form  (EPFBE).  The  monoenergetic  form  of  this  equation  can  be  repre¬ 
sented  as 


-A-  V  Km  A-  v'flf.  Ai  +  GyVc?,*.! 

=  -  A  V  A)  *  l  r,  A; 

(1.9) 

where  all  previous  definitions  apply  and  the  following  are  presented 

Y  *  2 

(even-parity  fluence) 

(1.10) 

>S^  =  2  ( Scf-Zl -•$(?,-*> ) 

(odd-parity  source) 

(1.11) 

•S^  =  2  (  -5<?  A)  *  Sit.-A)) 

(even-parity  source) 

(1.12) 

The  two  operators  Ku  and  Gg  are  positive  definite  and  self-adjoint  (Ref. 

17)  and  can  be  represented  as 

Km  -f(AJ  =  (  ~f(sL)  ■ 

►  /  A'-A){(/UcLa.') 

tTI 

(1.13) 

fai  =  <yr(t)  ~f<Ai  - 

j"  <ySj  ( r,  a'  A  1  (a!  1  dk 

HV 

(1.14) 

where 


0;u(r(  a'A) 


oo 

■iL 


<¥  ~Crl 


cr(?i  -  i 


C£u(rj  =  odd  Legendre  coefficients 
jP(A/U  =  /th  order  Legendre  polynomial 


Cl. 15) 


(1.16) 


O^A'AM  i  (q^<rA'/vJ  +  qi(^A'-Aj  ) 

(even-parity  differential  scattering  cross  section) 

The  vacuum  boundary  condition  becomes 

-  A  V  'f(r,A.)  )  =  O 

for  r  on  the  vacuum  surface  and  A- n  <  O 
For  the  reflective  boundary  condition 
'Pi r ,  a)  -  Ylr,A.-} 

where  r  is  on  the  reflective  boundary  and  Jd  is  defined  by  k-A-^k- a 
and CAx a'1 -n.  =  O  .  The  derivation  of  Eq.  (1.9),  Eq.  (1.13)  and  the 
boundary  condition  is  presented  in  Appendices  A  and  B  respectively. 

In  conjunction  with  this  self-adjoint  form,  a  functional  was 
derived  whose  Euler  equation  is  the  even-parity  fluence  equation  (Ref. 
17).  Kaplan  and  Davis  (Ref.  18)  used  this  functional  to  solve  a  simple 
monoenergetic  Milne  problem  with  isotropic  scatter.  They  derived  a 
coupled  set  of  differential  equations  that  represented  the  Euler  equa¬ 
tions  of  the  functional  for  the  various  regions  of  the  spatial  domain. 
This  set  of  equations  was  solved  by  finite  difference  methods.  The 
angular  dependence  of  the  neutron  distribution  in  their  problem  was 
represented  by  trial  functions  obtained  from  an  intuitive  knowledge 
of  the  true  solution.  The  use  of  such  predetermined  functions  is 
known  as  flux  synthesis  (Refs.  19-22)  and  provides  a  method  of  intro¬ 
ducing  a  priori  knowledge  of  the  particle  distribution  into  the  trial 
solution.  This  a  priori  knowledge  may  come  from  experience,  subsidiary 
calculations,  or  intuition. 
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Further  success  of  the  flux  synthesis  technique  was  achieved  by 
Roberds  and  Bridgman  (Ref.  23)  who  applied  the  flux  synthesis  method 
to  an  air-over-ground  problem  using. a.- weighted  residual  formulation  of 
the  Boltzmann  transport  equation.  Their  results  were  within  ten  percent 
of  those  generated  by  discrete  ordinates  and  reduced  the  computational 
time  by  a  factor  of  seven.  Also  significant  is  the  fact  that  their 
results  had  no  ray  effects,  since  the  chosen  synthesis  functions  changed 
the  angular  mesh  at  each  spatial  point  and  thus  did  not  allow  any  specific 
rays  to  penetrate  the  entire  spatial  domain. 

Miller,  Lewis,  and  Rossow  (Refs.  24  and  25)  applied  the  even- 
parity  fluence  functional  with  a  finite  element  trial  solution  in  space 
and  angle  to  solve  both  a  one  and  two-dimensional  monoenergetic  neutron 
transport  problem  in  plane  geometry.  In  the  2-D  problem  isotropic 
sources  and  scattering  were  used.  From  these  solutions,  the  authors 
concluded  that  the  finite  element  method  could  produce  results  with  a 
computational  time  savings  that  were  comparably  accurate  to  discrete 
ordinates  results.  Also,  as  in  the  case  of  synthesis  functions, 
finite  elements  used  with  the  EPFBE  are  not  subject  to  ray  effects.  A 
finite  element  method  has  been  incorporated  into  a  production  computer 
code  called  TRIPLET  (Ref.  26).  This  code  uses  the  finite  element  method 
to  represent  only  the  spatial  variation  of  the  neutron  distribution 
in  the  first-order  Boltzmann  equation.  The  angular  variation  is  repre¬ 
sented  by  discrete  ordinates  and  thus  this  code  is  subject  to  ray 
effects. 

In  1974  Kaper,  Leaf,  and  Lindeman  published  a  study  (Ref.  27)  which 
evaluates  several  finite  element  methods  for  solving  the  multigroup 
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neutron  transport  equation  using  the  variational  form  of  the  even-parity 
equation.  Their  final,  conclusions  were  unfavorable  and  lead  to  the 
recommendation  that  "the  use  of  high  order  approximation  procedures, 
based  on  finite  element  methods  and  applied  to  the  self-adjoint  form  of 
the  transport  equation  does  not  provide  a  viable  alternative  to  the  use 
of  the  discrete  ordinates  methods,  applied  in  combination  with  either 
finite  differences  or  finite  element  methods  to  the  standard  form  of  the 
transport  equation." 

These  authors  made  this  conclusion  based  on  their  experience  in  solv¬ 
ing  reactor  type  problems.  In  this  problem  ray  effects  do  not  occur  and 
thus  would  not  have  influenced  the  quoted  conclusion.  In  the  air-over¬ 
ground  problem  the  ray  effect  is  predominant  and  the  elimination  of  this 
effect  from  a  transport  solution  is  important.  Though  the  EPFBE  is  appli¬ 
cable  for  use  in  both  types  of  transport  problems,  its  value  in  mitigating 
or  eliminating  the  ray  effect  is  most  pronounced  in  an  air-over-ground 
problem.  Kaper,  Leaf,  and  Lindeman's  conclusion  does  not  apply  to  trans¬ 
port  problems  that  exhibit  the  ray  effect. 

In  1975  Briggs,  Miller,  and  Lewis  (Ref.  28)  did  an  extensive  theoret¬ 
ical  study  to  determine  the  reasons  behind  the  ray  effect  eliminating 
property  of  the  EPFBE.  Their  approach  was  to  solve  the  self-adjoint 
form  of  the  Boltzmann  equation  cast  in  a  variational  formulation  using 
three  different  angular  representations.  One  angular  formulation  was 
the  standard  discrete  ordinates  treatment  of  the  angular  domain.  The 
other  two  formulations  incorporated  either  piecewise  constant  or  piece- 
wise  bilinear  finite  elements  as  angular  trial  functions.  Solutions 
obtained  from  these  various  formulations  demonstrated  that  both  finite 
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element  angular  representations  eliminated  ray  effects,  while  the  dis¬ 
crete  ordinates  treatment  of  the  angular  variable  did  not.  To  explain 
their  computational  results  Briggs, -et  al ,  demonstrated  that  the 
operator  of  the  self-adjoint  form  of  the  Boltzmann  equation  changed 
from  a  hyperbolic  to  an  elliptic  form  when  the  finite  element  angular 
representation  replaced  the  discrete  ordinates  treatment.  This  change 
was  synonymous  with  the  disappearance  of  the  characteristic  lines  along 
the  discrete  ordinates  directions  of  neutron  travel  and  was  equivalent 
to  introducing  into  the  streaming  operators  fictitious  derivatives 
normal  to  these  lines  (Ref.  29).  These  derivatives  result  from  the 
averaging  of  the  transport  operator  over  the  solid  angle  support  of  the 
finite  element  basis  functions.  This  averaging  allows  a  coupling  of 
angular  directions  not  present  in  the  discrete  ordinates  method  and 
eliminates  the  ray  effect  from  transport  solutions. 

Purpose  of  the  Research 

The  second  order  even-parity  form  of  the  Boltzmann  transport  equa¬ 
tion  offered  a  sound  theoretical  basis  for  eliminating  the  ray  effect. 
However,  a  method  based  on  this  equation  is  much  more  difficult  to  pro¬ 
gram  than  the  discrete  ordinates  method  and  requires  the  added  expense 
of  evaluating  many  integrals.  Previous  applications  of  this  equation 
had  been  limited  to  one  dimensional  problems  with  linear  anisotropic 
scatter  and  two-dimensional  plane  geometry  with  only  isotropic  scatter. 
Extending  the  use  of  this  equation  to  two-dimensional  cylindrical  geom¬ 
etry  with  anisotropic  scatter  represented  a  significant  departure  from 
past  work.  The  purpose  of  this  research  was  to  determine  the  feasibility 
of  eliminating  the  ray  effect  when  solving  the  air-over-ground  problem  by 
using  the  EPFBE. 
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In  Chapter  II  the  Galerkin  method  is  applied  to  the  weak  form  of 
the  EPFBF.  Both  a  finite  element  and  synthesis  trial  solution  were  used. 
The  finite  element  trial  solution  was  a  tensor  product  consisting  of 
bilinear  Lagrange  polynomials.  The  synthesis  trial  solution  used  bi¬ 
linear  Lagrange  polynomials  to  represent  the  spatial  variation  of  the 
neutron  fluence  and  a  ellipsoidal  synthesis  function  expansion  to 
approximate  the  angular  dependence.  The  results  of  this  chapter  demon¬ 
strate  that  neither  the  Galerkin  or  Ritz  method  would  provide  a  feasible 
means  of  solving  the  air-over-ground  problem  since  the  resulting  infe- 
grals  cannot  be  efficiently  evaluated  either  numerically  or  analytically. 

In  Chapter  III  the  col  location  method  is  applied  to  the  EPFBE.  A 
trial  solution  consisting  of  cubic  splines  and  specially  derived  angular 
synthesis  functions  were  used.  This  method  allowed  the  analytic  integra¬ 
tion  of  the  resulting  integrals,  but  introduced  a  numerical  deficiency 
associated  with  the  odd-parity  fluence  transformation.  This  deficiency 
is  related  to  the  requirement  that  the  derivatives  with  respect  to  all 
phase  space  variables  of  the  selected  trial  solution  contain  the  neces¬ 
sary  information  to  accurately  represent  the  anisotropic  portion  of  the 
Boltzmann  fluence. 

In  Chapter  IV  the  problems  that  result  from  applying  a  multigroup 
method  to  the  EPFBE  are  outlined.  These  problems  center  on  the  odd-parity 
source  term  and  the  dependence  of  this  term  on  the  odd-parity  fluence. 

The  nested  integrals  that  result  are  complex  and  their  evaluation  would 
lead  to  problems  similar  to  those  encountered  in  Chapter  II. 

Chapter  V  presents  recommendations  and  conclusions. 


II.  Weak  Form  Equation  and  Galerkin  Method 


Introduction 


In  this  chapter  the  weak  form  of  the  EPFBE  is  formulated  for  the  air- 
over-ground  problem  using  the  Galerkin  method.  Evaluation  of  the  integrals 
that  result  from  this  method  is  attempted  using  both  analytical  and  numer¬ 
ical  integration  techniques.  Trial  solutions  consisting  of  linear  Lagrange 
polynomials  in  space  and  a  specially  derived  synthesis  function  in  angle 
are  used  to  approximate  the  even-parity  fluence.  Different  combinations 
of  the  above  integration  techniques  and  trial  solutions  are  attempted. 

The  results  demonstrate  that  the  Galerkin  and  Ritz  methods  introduce  too 
much  complexity  to  efficiently  assemble  a  global  matrix  representation 
of  the  EPFBE  for  an  air-over-ground  problem.  This  statement  applies  in 
general  to  the  Lagrange  polynomial  trial  solution  and  specifically  to 
the  synthesis  trial  solution. 

Weak  Form  and  Functional  Representations 

The  Galerkin  and  Ritz  methods  are  commonly  used  in  the  finite  ele¬ 
ment  method.  The  Galerkin  method  is  applied  to  the  weak  form  of  an  equa¬ 
tion.  The  weak  form  of  the  EPFBE  can  be  represented  as  (see  Appendix  C) 

X  (U-AVCfrTAJ,  (A  vYlt/i/)>+  <  €(K  A.I  >3  d?  = 

where  is  referred  to  as  a  weight  function.  V  represents  the 

spatial  domain,  Sv  is  the  part  of  the  surface  bounding  the  spatial  domain 
on  which  the  vacuum  boundary  condition  is  applied,  and 

<  CUr  A),  =  /  a(?.A>  bcAAjc/si  (2.2) 

L.  _ 
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The  Ritz  method  is  applied  to  a  functional  that  represents  an  equation. 


For  the  EPFBE  a  functional  representation  exists  and  can  be  written  as 

F  [Yc?,fil]  ~  (  a-jttl  £aV  Ycf'A.))  *  <TTlP)  '  (2 .3) 

/  .  'td 

*  [avYsaJ  a-  v  f (?,*.')  ct/i'  "  2  YifA.jSsl?,Aj 

-  4(7i  SucSJL)[a  ■  v  Yc?A  J  v  ¥( ?*ijt crsJt.A'AjSM(>:X)ok‘>v 


+  z«  V<.r,A)»> 


where 


<  aCr,sL)>v~  [  oL ?  )  uCf,A)  dA 

Y  'v  ■'Wn 

(2.4) 

<<  ait, A)  yy-  fdsfjA  n)  cU^A)dA 

(2.5) 

The  solutions  to  Eqs.  (2.1)  and  (2.3)  naturally  satisfy  the  vacuum  boundary 
condition.  The  reflective  boundary  condition  is  essential  for  both  equa¬ 
tions  and  must  be  imposed. 

If  the  same  trial  solution  is  used  in  solving  the  EPFBE,  the  Ritz 
and  Galerkin  methods  generate  an  identical  set  of  simultaneous  algebraic 
equations.  This  occurs  because  the  Ku  and  Gg  operators  in  this  equation 
are  both  positive  definite  and  self-adjoint.  The  Galerkin  method  was 
selected  for  use  in  this  chapter,  because  it  represented  a  convenient 
starting  point  for  determining  the  feasibility  of  applying  the  finite 
element  method  to  the  EPFBE. 

Galerkin  Method 

The  Galerkin  method  falls  under  the  broader  category  of  weighted 
residual  methods.  These  methods  assume  that  the  exact  solution  of  a 
differential  equation  is  not  known  and  must  be  approximated,  thus  intro¬ 
ducing  into  the  equation  an  error  called  the  residual.  The  weighted 
residual  method  requires  this  residual  to  vanish  in  some  average  sense 
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over  the  domain  of  definition  of  the  dependent  variable.  As  an  example 
consider  a  function  defined  over  some  domain  D  dependent  on  the 

variables  represented  by  the  vector  lx.. .  Let  the  behavior  of  <j>Cxl  in 
D  be  given  by 

L<j)(x)  =  fcx/  (2.6) 

where  -fcfJ  is  a  known  function  of  the  same  independent  variables.  Eq. 
(2.6)  can  be  rearranged  to  give 

L  d>c*;-  fat)  *  O  (2-?) 

Since  is  not  known,  an  approximate  solution  must  be  used.  Let  the 

approximate  solution  be  represented  by  .  Substituting  into 

Eq.  (2.7)  yields 

L  $(?)  -  f(?)  *  O  (2.8) 


or  stated  in  another  manner 

L  $OtJ  -  {(*)  =  R(X) 


(2.9) 


R(*J  is  called  the  residual  and  results  from  the  use  of  the  approximate 
solution  <$ct)  .  Multiplying  Eq.  (2.9)  by  the  weight  function  ccjt;  and 
integrating  over  the  domain  D  results  in 


f  [  f  cx;  ]  u dt  sfl Ret)  cci) dt 


(2.10) 


so 


Qti)  is  often  represented  by  a  linear  combination  of  a  linearly  inde¬ 
pendent  set  of  known  functions: 

M 

<pcx)  =  Z  OiQiC?;  (2.ii) 
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In  Eq.  (2.11)  the  a  's  are  the  constant  coefficients  and  the  Q  's  repre¬ 
sent  .  ,e  known  functions.  Eq.  (2.11)  is  called  a  trial  solution.  Insert 


ing  Eq  (2.11)  into  Eq.  (2.10)  gives 

}  [  L  Z  Q,  Ol  c s?>  -  -f( x>  J  ecx)  dx  -  j  Rci)  cc  x)  dx 


(2.12) 


Requiring  that  the  weighted  residual  vanish  (Ref.  30)  yields 
j  R(xJ  £(x)  dx  -  O- 


(2.13) 


Eq.  (2.13)  implies  that  the  weighted  residual  method  determines  a  set  of 
a  coefficients  that  requires  the  residual  function  to  be  orthogonal  to 
the  weight  function  and  thus  forces  the  residual  to  vanish  in  an  integral 
or  average  sense.  The  final  form  of  the  weighted  residual  method  is 

f  [  L  T.  CS)  -  -fcxiJ  £CxJ  d?  *  O 

J  n  i"  1 


If  the  weight  function  is  allowed  to  be  any  linear  combination  of  the 
same  known  function  as  the  trial  solution,  i.e. 

ecx)--  Z  LO;Cx)  (2.14) 

then  the  Galerkin  method  is  defined  and  requires 

f~  [  L  Z  Qi  O’  Cxi  -  {(ij  J  Qj(x)  dx  =  -  K?,  •  M  (2.15) 

Transport  problems  are  normally  of  such  complexity  that  an  analytical 
solution  is  not  possible.  These  problems  necessitate  the  use  of  a  trial 
solution  that  describes  both  the  spatial  and  angular  variation  of  the 
neutron  fluence.  Such  a  trial  solution  for  the  even-parity  fluence  can 
be  presented  as 
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(2.16) 


W  M 


-  Z  Z  ai;SLL^)  At  A) 


where  S±t*i  represents  a  known  function  to  approximate  the  spatial 
variation  and  Ajt*)  plays  a  similar  role  for  the  angular  dependence. 
If  the  Galerkin  method  is  used  with  the  trial  solution  in  Eq.  (2.16), 

then  the  proper  form  of  the  weight  function  is 

w  n 

1  Z  bKL  SKC?)  Al(a.) 

KS  I  t«  • 


(2.17) 


Substituting  Eqs.  (2.16)  and  (2.17)  into  the  weak  form  of  the  EPFBE, 
Eq.  (2.1),  yields 


ft 

«•»'  ia<  J  V 


<-A.  v  AL(X),  A  V  Stifi)  Aj  (A)  >  + 


(2.18) 


SKt?)  Alu)'  5t c?jAj(*)>)dt+  sK(t)AL(A)t 5V; Aj(a ) Ia  a bds) 


-ft 


SK(vAL(A)tKM$Mct;Ai>  +  <SKc*}ALiA]tS5c?,*/>]  d? 

’for’  1C  5  /,  2,  •  •  *  *  |  ** 

Ar  t  4  h  2,  .  ••  •  t  M 


Eq.  (2.18)  is  referred  to  as  the  Galerkin  weak  form  equation  throughout 
the  remainder  of  this  dissertation. 

Numerical  Integration 

The  trial  solution  initially  used  with  Eq.  (2.18)  consisted  of 
piecewise  bilinear  Lagrange  polynomials  to  approximate  both  the  spatial 
and  angular  variation  of  the  fluence.  These  polynomials  have  only  C 

o 

continuity  (continuity  of  the  function  only).  This  degree  of  continuity 
is  all  that  is  required  since  the  highest  derivative  in  the  Galerkin 
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weak  form  equation  equals  one  (Ref  31).  In  any  given  element,  the  form 
of  these  polynomials  is  given  by 

L;  Cx.yj  =  Q;  X y  +  &-X.  +•  Q y  ♦  oil  (2*19) 

where  the  i  subscript  refers  to  a  particular  node  in  the  discretized 

domain  consisting  of  rectangular  elements.  Since  four  nodal  degrees  of 

freedom  exist  in  each  rectangular  element,  the  coefficients  in  Eq.  (2.19) 

can  be  determined.  These  coefficients  satisfy  the  requirements  that  at 

* 

one  node  L^C*. y;  equals  one  and  at  the  other  three  the  polynomial  is 
zero.  These  constraints  lead  to  the  classic  tent  function  illustrated 
in  Fig.  3.  The  same  procedure  is  followed  in  determining  the  other  tent 
functions  for  the  remaining  nodes  in  the  rectangular  element  and  also 
the  remaining  elements  in  the  mesh.  This  procedure  is  given  in  more 
detail  in  Appendix  D. 

The  Lagrange  polynomials  representing  the  spatial  and  angular  varia¬ 
tion  of  the  fluence  are  combined  using  a  tensor  product.  In  this  tensor 
product,  the  spatial  and  angular  dependence  of  the  fluence  is  assumed 
separable  in  a  manner  similar  to  the  Pn  method  (Ref.  1).  Using  the  Galerkin 
method  with  this  bilinear  Lagrange  polynomial  tenser  product  trial  solution, 
one  defines  the  following  representations  for  the  even-parity  fluence  and 
weight  function  respectively: 


.  H  N 

Ur.Aj  3  2  Z  CU:L<?/  L:UJ 

(2.20) 

€Cr,A)-  Z  t  bKL  Lk(*jLl(Z) 

IC«I 

(2.21) 

where  M  represents  the  total  number  of  spatial  nodes  and  N  the  number 
of  angular  nodes. 
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The  trial  solution  presented  in  Eq.  (2.20)  was  used  (Ref.  26)  to 
find  the  solution  of  a  monoenergetic  transport  problem  involving  two-dimen¬ 
sional  plane  geometry  and  isotropic  scatter.  The  success  of  this  appli¬ 
cation  depended  in  part  on  the  simplicity  that  the  isotropic  scatter 
assumption  introduced.  Scattering  from  one  angular  direction  to  another 
is  represented  in  the  first  order  Boltzmann  equation  by  the  differential 
scattering  cross  section  oj(r; a'-*-* A)  for  scattering  from  A  to  A  . 

The  differential  scattering  cross  section  is  normally  approximated  by  an 
expansion  in  terms  of  Legendre  polynomials,  i.e. 

p 

crs  ( r  i  (2.22) 


where  is  the  JC"  order  Legendre  polynomial,  is  a  spatially 

dependent  coefficient,  and  p  represents  the  order  of  the  approximation. 
In  the  case  of  isotropic  media  the  A  -*/ 1  variable  is  replaced  with 
which  represents  the  cosine  of  the  angle  between  A'  and  A  and  is 


expressed  as 


/JL4  * 


uw.'  + 


Co.5  (X-  X' ) 


(2.23) 


For  the  EPFBE  the  definition  presented  in  Eq.  (1.15)  and  (1.16)  leads  to 

p 


03<}(r,Uf)  -  ^  OsjilrJ  p£(.Uj) 


(2.24) 


2  = 


(2.25) 


The  assumption  of  isotropic  scattering  reduces  the  complexity  of  Eq.  (2.24) 
and  (2.25)  to 

*  Qi#  tr  ;  (2.26) 
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(2.27) 


In  an  air-over-ground  problem,  the  differential  scattering  cross  sec¬ 
tion  is  expressed  as  a  third  order  Legendre  expansion  (P=3).  The 
introduction  of  anisotropic  scattering  significantly  complicates  Eq.  (2.18) 
and  the  process  of  reducing  this  equation  to  a  matrix  form.  The  complica¬ 
tion  that  occurs  in  the  Gg  and  Ku  operators  can  be  seen  by  comr_ring  the 
isotropic  scattering  case  with  an  anisotropic  representation.  Combjjiing 
the  definition  of  these  operators  presented  in  Eqs.  (1.13)  and  (1.14)  with 
the  isotropic  scatter  results  of  Eqs.  (2.26)  end  (2.27),  leads  to  the 
following  representation  of  these  operators  in  Eq.  (2.1): 


Tr(.n'f(.rA)  -  < 7s4Crjl  Yc 


(2.28) 


1 

aft? 1 


[ 


*  — * 
A*  7 


(2.29) 


Eqs.  (2.28)  and  (2.29)  demonstrate  two  simplifications  that  result  from 
the  isotropic  scatter  assumption.  The  first  is  the  absense  of  the 
integral  term  in  the  Ku  operators,  since  *  °  .  The  second  i  s 

the  simple  form  of  the  integral  in  the  Gg  term  that  results  from  the  P0 
Legendre  polynomial  equaling  one.  Using  the  trial  solution  of  Eq.  (2.20), 
this  global  integral  becomes  an  evaluation  of  the  bilinear  Lagrange  poly¬ 
nomials  over  the  angular  domain.  Since  these  polynomials  are  identically 
defined  over  each  angular  element,  the  integral  can  be  evaluated  using 
the  results  of  a  single  canonical  element  evaluation.  The  value  of  this 
integral  for  any  angular  element  can  be  found  by  simply  multiplying  the 
area  of  the  angular  element  by  the  canonical  result.  (See  Appendix  E 
for  more  detail.)  If  a  more  precise  value  is  necessary,  only  the  single 
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canonical  evaluation  must  be  repeated.  This  simple  procedure  Is  possible 
due  to  the  definition. of  this  integral. 

The  integrals  resulting  from  both. the  Gg  and  Ku  operators  couple  all 
scattering  directions  In  the  angular  domain  and  thus  must  be  evaluated 
over  all  angular  finite  elements.  This  global  definition  is  contrary 
to  the  other  integrals  resulting  from  the  definition  of  the  weak  form 
equation.  These  integrals  are  defined  locally  over  a  single  spatial 
and  angular  finite  element  and  during  the  assembly  of  the  global  matrix 

m  * 

are  evaluated  over  a  single  element  at  a  time.  The  difference  between 
a  local  and  global  definition  can  be  clarified  by  examining  the  term  in 


Eq  (2.1)  that  contains  the  Gg  operator. 


j'  EC^AJ  cAA' J ci?  (2.30) 


To  determine  the  value  of  one  element  in  the  global  matrix,  the  two  outer 
integrals  of  Eq.  (2.30)  are  evaluated  over  a  few  spatial  and  angular 
finite  elements.  For  the  same  single  value  the  inner  integral  must  be 
evaluated  over  the  entire  angular  domain.  Since  the  neutron  fluence  is 
highly  anisotropic  near  the  source  in  an  air-over-ground  problem  the 
angular  domain  requires  a  refined  discretization  to  accurately  approximate 
this  type  of  distribution.  Thus  evaluating  the  global  integrals  result¬ 
ing  from  the  and  Gg  operator  could  be  quite  costly,  if  the  definition 
of  the  integrand  prohibited  the  effective  use  of  canonical  integration. 
Such  is  the  case,  when  an  anisotropic  scattering  process  is  used. 


The  Gg  and  Kg  operators  take  within  Eq.  (2.1)  the  following  forms 
in  the  anisotropic  case: 


G3  fi^A)  ~  crT(r)  YirAl-fcqir.A  A)  Vlr.A'jclA’  (2.31) 
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KM/V  ?  Wf-Aj  =  a,*,  ( A vftt.XialA-  )  (2.32) 

JHn 


The  degree  by  which  these  forms  deviate  from  the  isotropic  case  depends 
on  the  order  of  the  even  and  odd-parity  differential  scattering  cross 
sections.  To  illustrate  the  complications  introduced  by  anisotropic 
scattering  in  these  terms,  assume  a  first  order  scattering  representation 
involving  the  PQ  and  P-j  terms  of  the  Legendre  expansion.  This  order  of 
scattering  approximation  does  not  alter  the  definition  of  the  Gg  teem 
presented  in  Eq.  (2.28),  but  does  add  the  integral  term  to  the  Ku  opera¬ 
tor,  i.e. 


Substituting  the  definition  of  into  Eq.  (2.33)  results  in  the  follow¬ 
ing  formulation  for  the  global  integral  term. 

Xif,*.)*  CS.tfJ  (£  i  dx'Uu.'  t  (2.34) 

\fi  ~ ^ 1 'coil*- x  iA'  vYi t\‘)  dt'cljj.') 

The  two  integrals  in  Eq.  (2.34)  are  considerably  more  complex  than  those 

resulting  from  the  Gg  operator  due  to  the  presence  of  the  gradient  term. 

The  increased  complexity  of  these  integrands  affect  the  order  of  quadra¬ 
ture  needed  to  accurately  evaluate  them  and  thus  the  computational  time 
required.  Another  significant  deviation  from  the  isotropic  case  is  the 
presence  of  \J i -/*•*'  and  co*(x-x‘l  within  these  integrands.  These  func¬ 
tions  cause  the  integrand  to  be  defined  differently  in  each  angular 
finite  element,  thus  altering  the  basic  premise  that  allowed  the  efficient 
evaluation  of  the  integral  in  the  isotropic  case.  The  presence  of  these 
functions  in  the  integrand  now  force  the  global  integral  to  be  evaluated 
over  each  angular  finite  element,  rather  than  a  single  canonical  element. 
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Since  this  global  evaluation  must  be  repeated  for  each  angular  quadrature 
point,  the  expense  of  .forming  a  global  matrix  representing  the  Galerkin 
weak  form  equation  is  greatly  increased. 

The  problems  demonstrated  with  a  P-j  anisotropic  scattering  representa¬ 
tion  for  the  «u  operator  become  more  severe  with  a  P3  representation  and 
spread  to  the  global  integral  term  of  the  Gg  operator.  This  Gg  operator 
is  affected  by  the  ?2  term  Legendre  expansion  that  introduces 

\JT1/Z'77'  and  CosU-%'  )to  the  integrand  of  this  integral.  The  Legendre 
expansion  of  the  angular  scattering  significantly  increases  the  com¬ 
plexity  of  an  anisotropic  scattering  problem  compared  to  an  isotropic 
scattering  problem.  This  complexity  is  primarily  concentrated  in  the 
global  integrals  resulting  from  the  Gg  and  Ku  operators. 

A  modified  computer  code  (Ref.  32)  was  used  to  determine  the  exact 
effect  anisotropic  scattering  had  on  computer  execution  time  and  accuracy. 
This  code  assembled  the  global  matrix  resulting  from  the  Galerkin  weak 
form  equation  Eq.  (2.18)  using  a  bilinear  Lagrange  tensor  product  trial 
solution.  The  evaluation  of  all  integrals  resulting  from  this  equation 
was  performed  during  assembly  using  numerical  integration.  This  integra¬ 
tion  was  carried  out  over  a  canonical  representation  of  the  spatial  and 
angular  elements.  A  flow  diagram  of  this  code  appears  in  Appendix  H. 

The  problem  selected  to  examine  the  effect  of  anisotropic  scattering 
was  a  simple  air  burst  problem.  The  spatial  domain  consisted  only  of  air 
and  was  discretized  into  a  single  rectangular  finite  element  extending 
from  0  to  100  meters  in  both  the  £>  and  l  directions.  The  angular  domain 
was  zoned  into  a  single  rectangular  finite  element  with  the  u  variable 
extending  from  0  to  1  and  X  from  0  to7T.  The  symmetry  of  the  even-parity 
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fluence  allowed  this  reduced  angular  domain.  All  cross  sections  were 
taken  from  the  first  energy  group  of  . the  DLC-31  set  generated  by  the 
Radiation  Shielding  Information  Center,  Oak  Ridge,  Tennessee  (Ref.  34). 

The  vacuum  boundary  condition  was  satisfied  naturally  on  three  boundaries 
and  the  reflective  boundary  condition  was  enforced  on  the  cylinder  axis. 

The  global  matrix  resulting  from  the  definition  of  this  problem  was 
16  x  16. 

Four  problems  were  run  on  the  CDC-6600  computer  using  this  code. 

These  problems  differed  only  in  the  order  of  the  scattering  approximation 
and  quadrature  set  used.  The  effect  of  increasing  the  order  of  the  scat¬ 
tering  representation  was  measured  by  two  parameters.  One  parameter  was 
the  computational  time  required  to  compute  and  assemble  a  global  matrix. 

The  second  parameter  was  concerned  with  the  sign  of  the  eigenvalues  related 
to  this  matrix.  The  sign  of  all  eigenvalues  of  a  global  matrix  from  the 
Galerkin  weak  form  equation  should  be  positive.  This  results  from  the 
positive  definiteness  of  the  and  operators  used  in  this  equation. 

A  negative  eigenvalue  would  indicate  that  the  elements  contained  within 
a  global  matrix  were  inaccurately  computed.  Since  all  computations  in 
the  Galerkin  weak  form  equation  center  around  the  evaluation  of  integrals, 
this  would  indicate  that  the  order  of  the  quadrature  set  used  was  inadequate 
and  would  need  to  be  Increased.  Increasing  the  order  of  the  quadrature 
would  increase  the  computational  time  required  to  assemble  a  global  matrix. 
Thus,  it  can  be  seen  that  the  two  evaluation  parameters  are  related  and 
together  indicate  the  relative  efficiency  of  assembling  an  accurate  global 
matrix  representation  of  the  Galerkin  weak  form  equation. 


The  results  of  the  four  sample  problems  are  presented  in  Fig.  4  and 
substantiate  the  analysis  previously  presented.  Problem  one  is  analogous 
to  the  analysis  that  used  an  isotropic  .scattering  representation.  Here 
a  simple  two  point  Gaussian  quadrature  set  was  adequate  to  generate  an 
accurate  global  matrix.  Problem  two  and  three  demonstrated  the  increasing 
complexity  that  occurred  as  the  order  of  the  scatter  approximation  increased. 
For  both  problems  a  two  point  Gaussian  quadrature  set  was  used.  The 

increased  complexity  was  reflected  in  the  number  of  negative  eigenvalues 

'  , 

associated  with  the  respective  global  matrices.  Problem  four  was  an 
attempt  to  increase  the  order  of  the  quadrature  set  for  a  P=1  scattering 
approximation.  This  problem  demonstrated  that  changing  the  order  of  the 
quadrature  set  substantially  increased  the  computational  time  required 
to  generate  a  global  matrix.  This  occurred  because  of  the  six  nested 
integrals  that  are  contained  in  the  Galerkin  weak  form  equation.  If 
Eq.  (2.33)  is  substituted  into  the  first  term  of  Eq.  (2.1)  the  following 
would  result 


Ji  J  j  -A.’  ^  €(r.  A.)  ( /\.  V  Vc^A.1  *■  'll  ) clxcli  oL/a  cLp  (2.35) 

C  V-  •'fc  * 


Eq.  (2.35)  clearly  illustrates  these  nested  integrals.  Evaluating  these 
integrals  requires  Q®  evaluations  per  matrix  element,  where  Q  represents 
the  number  of  quadrature  points  for  each  variable.  For  a  two  point  quadra¬ 
ture  set  64  evaluations  are  needed  per  matrix  element  and  four  quadrature 
points  would  require  4096  evaluations  per  element.  Considering  the  256 
elements  in  the  matrix  the  substantial  increase  in  computational  time  is 
not  surprising. 

The  numerical  integration  procedure  used  in  the  computer  code  that 
generated  the  results  presented  in  Fig.  4  was  not  efficient.  An  analysis 
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Fig  4.  Results  Using  Galerkin  Weak  Form  Equation 
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of  the  second  term  in  Eq.  (2.35)  demonstrated  it  could  be  separated  into 
9  spatial  and  27  angular  integrals.  The  evaluation  of  these  integrals 
using  a  four  point  quadrature  set  would  require  576  evaluations  per  matrix 
element.  Though  the  576  evaluations  represented  a  significant  reduction 
from  the  4096  previously  required,  it  still  was  not  enough  reduction  to 
support  the  use  of  numerical  integration. 

Certain  requirements  of  the  air-over-ground  problem  oppose  the  cost 
efficient  implementation  of  a  finite  element  solution  technique  usigg 
numerical  integration.  This  transport  problem  requires  the  use  of  a 
third  order  Legendre  expansion  to  approximate  the  scattering  process. 

This  expansion  used  in  the  weak  form  of  the  EPFBE  would  significantly 
increase  the  complexity  of  the  resulting  integrals.  As  demonstrated  in 
Fig.  4,  complex  integrals  require  the  use  of  higher  order  quadrature  sets 
to  assure  the  positive  definiteness  of  a  global  matrix.  The  air-over¬ 
ground  problem  also  requires  a  very  finely  zoned  angular  and  spatial 
domain.  Such  zoning  intensifies  the  problems  associated  with  evaluating 
the  globally  defined  angular  integrals  and  increases  the  number  of  local 
element  integrations  that  must  be  done.  Since  the  angular  integrals 
resulting  from  the  weak  form  of  the  EPFBE  can  not  be  efficiently  evalu¬ 
ated  by  canonical  techniques,  a  computational  cost  inefficiency  results. 
The  combined  effect  of  high  order  quadrature  sets  and  refined  meshing 
significantly  offset  even  substantial  gains  made  in  improving  the 
efficiency  of  a  numerical  integration  process.  For  this  reason  further 
work  based  on  a  finite  element  method  using  numerical  integration  was 
stopped  and  the  use  of  analytical  integration  techniques  were  attempted. 

The  success  of  previous  work  using  the  EPFBE  functional  was  based 
on  the  isotropic  scattering  assumption.  This  assumption  automatically 
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resolved  all  the  problems  that  rendered  either  this  formulation  or  the 
weak  form  equation  too  complex  for  efficient  numerical  solution.  - 


The  numerical  integration  problems  encountered  using  the  Galerkin 
weak  form  equation  are  not  unique.  Strang  and  Fix  (Ref.  33)  devote  an 
entire  section  in  a  chapter  entitled  "Variational  Crimes"  to  this  problem. 
In  this  section  they  theoretically  discuss  problems  similar  to  those 
presented  in  this  chapter  and  conclude  that  "it  is  very  important  to 
control  properly  the  fraction  of  computer  time  spent  on  numerical  integra- 
tion."  As  evidenced  by  the  sample  problems,  the  computer  time  necessary 
to  evaluate  the  integrals  resulting  from  the  Galerkin  weak  form  equation 
has  not  been  properly  controlled. 

Analytical  Integration 

A  natural  solution  to  the  numerical  integration  problem  was  to 
attempt  the  analytical  evaluation  of  the  integrals  associated  with  the 
Galerkin  weak  form  equation.  This  method  of  integration  would  alleviate 
both  the  accuracy  problem  and  computational  time  problems  met  in  the 
previous  section.  Analytical  integration  would  be  precise  and  require 
only  a  single  function  evaluation  per  finite  element,  instead  of  the 
many  evaluations  required  by  numerical  integration.  The  structure  of 
the  integrals  resulting  from  this  equation  were  known  to  be  complicated; 
however,  it  was  horcJ  that  the  MACSYMA  (Ref.  34)  symbolic  algebraic  mani¬ 
pulating  system  would  provide  the  necessary  simplification.  Two  differ¬ 
ent  trial  solutions  were  used  in  attempting  this  analytical  integration. 

The  first  trial  solution  was  the  tensor  product  consisting  of 
bilinear  Lagrange  polynomials  used  in  the  previous  section.  The  first 
term  on  the  right-hand  side  of  Eq.  (2.18)  was  input  to  the  MACSYMA  pro¬ 
gram  by  defining  the  individual  terms  that  composed  it.  As  an  example. 


the  /U*  variable  represented  by  Eq.  ( 2.23)  was  defined  in  MACSYMA.  Next 
the  and  Legendre  polynomials  were  defined  in  terms  of  yu^  .  Using 
the  preceding  two  definitions  the  odd^parity  differential  scattering 
cross  section  was  constructed.  The  Lagrange  polynomials  representing 
both  the  spatial  and  angular  variation  of  the  neutron  fluence  were  next 
defined.  In  the  first  term  of  the  Galerkin  weak  form  equation  both  the 
weight  function  and  trial  solution  is  operated  upon  by/v'b  .  In  cylin¬ 


drical  geometry  this  operator  is  defined  as 

CosCt I  f(*Xl  ) 


MACSYMA  was  programmed  to  perform  the  operation  defined  in  Eq.  (2.36)  on 
both  the  previously  defined  weight  function  and  trial  solution.  With 
all  the  individual  components  defined,  MACSYMA  assembled  the  first  term 
of  the  Galerkin  weak  form  equation. 

Several  different  techniques  were  attempted  to  perform  analytical 
integration  over  the  required  six  nested  integrals.  One  attempt  involved 
integrating  the  whole  assembled  term.  This  attempt  was  only  successful  in 
performing  the  two  intermost  integrations  over  the  u'  and  x'  variables. 
Further  integration  was  impossible,  since  the  number  of  terms  generated 
by  the  successful  integrations  exceeded  the  memory  capacity  of  MACSYMA. 
This  large  number  of  terms  developed  because  of  the  indefinite  integration 
limits  and  the  complex  nature  of  the  trigonometric  functions  that  result 
from  the  assembled  term.  The  integration  limits  were  arbitrarily  set 
so  the  dimensions  of  the  finite  elements  in  the  angular  and  spatial  domain 
could  be  varied  without  repeating  the  work  on  MACSYMA.  The  use  of  these 
arbitrary  limits  always  generated  two  terms  for  each  integration,  since 
an  upper  and  lower  limit  evaluation  was  necessary.  The  complex  trigono¬ 
metric  terms  resulted  from  the  product  of  those  involved  in  Eq.  (2.36) 
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and  the  definition  of  .  Integrals  involving  these  terms  required 
reduction  techniques  for  evaluation.  These  reduction  techniques  always 
generated  several  additional  terms.  The  exact  number  of  the  terms  was 
dependent  on  the  power  of  the  exponent  associated  with  the  trigonometric 
function. 

To  resolve  the  memory  capacity  problem,  several  different  MACSYMA 
functions  that  optimized  the  use  of  memory  were  tried.  The  use  of  these 
functions  did  not  solve  the  memory  problem  thus  motivating  the  use  of  a 
new  approach. 

The  previous  approach  attempted  to  evaluate  the  entire  first  term 
of  the  Galerkin  weak  form  equation,  while  this  new  approach  sought  to 
simplify  the  expression  into  a  number  of  smaller  less  complex  terms.  To 
accomplish  this  goal  a  special  MACSYMA  function  was  written  that  auto¬ 
matically  searched  the  string  of  simplified  terms  to  identify  those  that 
have  the  same  form  of  variable  dependence.  Once  identified  the  coeffi¬ 
cients  of  these  terms  were  added  together,  thus  reducing  the  total  number 
of  terms  requiring  evaluation.  Several  attempts  were  made  to  implement 
this  approach  and  use  this  specially  derived  MACSYMA  function.  In  all 
cases  memory  capacity  again  was  exceeded  during  the  simplification  process 
A  careful  examination  of  the  first  term  of  the  Galerkin  weak  form  equation 
demonstrates  this  expression  can  be  simplified  into  over  29,000  distinct 
terms  (see  Appendix  F).  MACSYMA  does  not  have  the  capacity  to  generate 
this  number  of  terms  and  no  successful  means  was  found  to  generate  only 
a  small  portion  of  these  terms  at  one  time. 

These  two  approaches  failed  for  reasons  similar  to  those  found  in 
attempting  numerical  integration.  Once  again  the  complexity  of  the  inte¬ 
grands  contributed  to  the  failure  of  using  the  Galerkin  weak  form  equation 


The  integral  evaluation  of  these  functions  necessiated  the  use  of  reduction 
formulas  that  required  the  use  of  significant  amounts  of  memory  in  MACSYMA. 
The  nested  integral  problem  also  surfaced.  In  the  analytical  integration 
attempt  this  problem  manifested  itself  by  taxing  the  limited  memory 
resources  available  to  MACSYMA.  Therefore,  any  means  that  would  reduce 
the  requirement  for  this  memory  would  improve  the  chances  of  successfully 
using  this  program. 

The  space  angle  synthesis  ( SAS )  method  discussed  in  Chapter  I  offered 
the  possibility  of  reducing  these  memory  requirements.  The  Lagrange  poly¬ 
nomial  trial  solution  used  in  an  air-over-ground  problem  requires  sixteen 
distinct  terms.  Substituting  a  global  synthesis  function  representation 
for  the  angularly  dependent  Lagrange  polynomials  would  condense  the  form 
of  the  trial  solution.  This  reduction  would  result  from  replacing  the 
four  terms  of  the  bilinear  Lagrange  representation  with  a  single  term. 

The  choice  of  synthesis  functions  was  influenced  by  the  work  of  Roberds 
and  Bridgman  (Ref.  23).  These  authors  used  an  off-centered  ellipsoidal 
function  to  represent  the  angular  variation  of  the  neutron  fluence  in 
an  air-over-ground  problem.  The  most  significant  feature  of  this 
ellipsoidal  function  was  its  ability  to  be  analytically  integrated  in 
a  weighted  residual  formulation  of  the  first  order  Boltzmann  equation. 
Additionally,  the  effectiveness  of  this  function  was  demonstrated  by 
accurately  representing  the  angular  variation  of  the  neutron  fluence  and 
generating  ray  effect  free  scalar  fluence  solutions. 

The  off-centered  ellipsoidal  function  used  by  Roberds  and  Bridgman 
needed  modification  for  use  in  the  Galerkin  weak  form  equation.  The 
even-parity  fluence  as  defined  in  Eq.  (1.10)  is  an  even  function  in  the 
angular  variable.  This  property  forces  the  ellipsoid  to  be  centered. 
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The  expression  for  a  centered  ellipsoid,  symmetric  about  the  position 


vector  from  the  point  source  to  position  )  Is  (see  Fig.  5) 

El£-A^  =  £>  V  I-  f (/*,*)' 


(2.37) 


where 


E(A  )  is  the  distance  from  to  the  ellipsoidal  surface, 

b  denotes  the  length  of  the  minor  axis,  tcpni  =  the  cosine  of  the  angle 
between  R  and  the  axis  of  revolution 

fc^.  v)  =  Coj TT-fjtT  cosixl 
The  r\_  parameter  is  now  a  function  of  the  (f.*  )  coordinates,  since  it 
is  dependent  on  the  spatial  location.  This  dependence  is  defined  in 
the  following  manner: 

cp- 

SlH(fl)  S  is)1’ 


( 1  -  ) _ 

C05''1’  '■  f.-v* 


where  the  coordinates  )  represent  the  spatial  location  of  the 

point  source  of  weapon  radiation.  In  the  air-over-ground  problem 
since  the  point  source  is  on  the  axis  of  the  cylinder. 

The  angularly  dependent  part  of  the  trial  solution  Eq.  (2.20)  was 
modified  to  include  these  ellipsoidal  synthesis  functions.  The  trial 
solution  chosen  replaced  the  angularly  defined  bilinear  Lagrange  poly¬ 
nomials  with  a  three  term  ellipsoidal  expansion 

M  3 


(2.38) 


where  3^)  is  still  the  spatially  defined  bilinear  Lagrange  polynomials. 
The  three  ellipsoidal  functions  were  chosen  to  represent  the  expected 
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spectrum  of  angular  variation  that  would  exist  for  an  air-over-ground 

problem  (see  Fig.  6). -  E,(A.I  was  chosen  to  represent  the  extremely 

anisotropic  nature  of  the  angular  fluence  near  the  point  source. 

represented  the  isotropic  nature  of  the  angular  fluence  at  large  distances 

from  the  source.  Ezia.)  was  chosen  to  act  as  a  transitional  shape 

between  the  two  others.  The  b  parameter  in  Eq.  (2.37)  determined  these 

shapes.  The  mixing  of  these  shapes  is  accomplished  by  the  coefficients. 

The  value  of  these  coefficients  is  determined  by  solving  the  olooal 

<*>  * 

matrix  representing  the  weak  form  of  the  even-parity  equation. 

Analytical  integration  of  these  ellipsoidal  synthesis  functions 
in  the  weak  form  of  the  even-parity  equation  was  not  possible.  This 


can  be  illustrated  by  examining  the  simplest  term  in  Eq.  (2.1). 

^  cLf* 

Substituting  Eq.  (2.38)  and  the  proper  definition  of  cctf/l;  for  the 

Galerkin  method  results  in: 
m3  r 

Z  Z  at; 

«.=  i  J*«  ’JJ?  K  L  '  J  <-  J 

for  K«  ■,•  •••  .  M 

l*  «  1,2,3 


(2.39) 


(2.40) 


If  isotropic  scatter  Is  assumed  the  following  angular  integral  would 
result  from  Eq.  (2.40) 

)  EL(vE:CA.tcM 

JHtt  J 


By  substituting  Eq.  (2.37)  into  the  above,  the  following  results: 


(2.41) 


36 


AMtSOTROWC 


Representative  Angular  Fluence  Shapes  (Polar  Plots) 


(2.42) 


Eq.  (2.42)  could  not  be  analytically  .integrated  by  the  MACSYMA  program. 

The  weighted  residual  technique  applied  to  the  first  order  Boltzmann 
equation  by  Roberds  and  Bridgman  never  produced  a  multiplication  of  two 
ellipsoidal  functions.  This  multiplication  occurred  because  the  Galerkin 
method  was  used. 

Numerical  integration  of  the  ellipsoidal  trial  function  was  also 
attempted.  The  computer  program  previously  described  was  modified  to 
replace  the  bilinear  Lagrange  polynomials  defined  over  the  angular 
domain  with  the  ellipsoidal  functions.  All  derivatives  of  this  function 
needed  for  an  evaluation  of  the  Galerkin  weak  form  equation  were  analyti¬ 
cally  calculated  and  formulated  in  FORTRAN  by  MACSYMA.  Sample  problems 
using  this  modified  code  demonstrated  the  need  for  high  order  quadrature 
sets  to  accurately  evaluate  the  resulting  integrals  and  thus  properly 
define  the  global  matrix.  As  before,  the  use  of  these  higher  order  quadra¬ 
ture  sets  made  the  computational  time  necessary  to  evaluate  one  spatial 
finite  element  prohibitive  for  an  air-over-ground  problem. 

Further  efforts  to  use  the  Galerkin  weak  form  equation  to  obtain 
solutions  to  an  air-over-ground  problem  were  not  attempted.  The  complexity 
of  this  equation  formulated  for  this  problem  had  been  demonstrated. 

Neither  numerical  nor  analytical  integration  could  efficiently  form  an 
accurately  defined  global  matrix  representation. 

The  next  logical  step  in  attempting  to  generate  a  solution  to  the 
EPFBE  was  to  use  a  less  complex  numerical  method.  The  numerical  method 
selected  was  collocation  as  described  in  the  next  chapter. 
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III.  Collocation  Method 


Criteria 

The  results  obtained  from  the  previous  chapter  identified  two  criteria 
that  must  be  met  to  successfully  generate  a  global  matrix  representation 
of  the  EPFBE  for  an  air-over-ground  problem.  The  first  criterion  is  to 
use  a  numerical  method  that  allows  the  efficient  computation  of  this 
matrix.  The  Galerkin  method  failed  this  criterion  by  requiring  the  use 
of  a  weight  function  that  significantly  added  to  the  complexity.  Addi- 
tionally  this  method  generated  the  nested  integrals  that  contributed 
to  the  problems  in  attempting  both  analytical  and  numerical  integration. 

The  second  criterion  is  to  use  a  trial  solution  that  is  simply  defined 
and  analytically  integrable.  The  bilinear  Lagrange  polynomials  were 
simple  in  form;  however,  they  expanded  the  number  of  distinct  terms  in 
the  Galerkin  weak  form  equation  beyond  the  memory  requirements  available 
in  MACSYMA.  The  ellipsoidal  angular  synthesis  function  had  only  one 
term  but  proved  not  analytically  integrable.  Also  the  complexity  of 
this  synthesis  function  prohibited  the  efficient  use  of  numerical 
integration,  since  high  order  quadrature  sets  were  needed  for  accurate 
evaluation. 

The  main  purposes  of  this  chapter  are  to  present  a  numerical  method 
that  satisfies  the  first  criterion,  to  develop  a  trial  solution  that 
meets  the  requirements  of  the  second  criterion,  and  to  identify  the 
problems  that  result  from  applying  this  combination  to  an  air-over-ground 
problem. 

The  numerical  technique  chosen  is  the  collocation  method.  This 
method  was  selected  after  reviewing  the  work  of  Houstis,  Lynch,  Rice, 
and  Papatheodorou  (Ref.  35).  These  authors  generated  solutions  to  seven¬ 
teen  different  elliptic  partial  differential  equations  using  the  numerical 
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methods  of  finite  differences,  collocation,  least  squares,  and  Galerkin. 

The  purpose  of  this  paper  was  to  compare  the  efficiency  and  accuracy  of 
these  numerical  methods.  This  research  indicated  that  the  collocation 
method  offered  an  efficient  and  accurate  alternative  to  the  other  three 
techniques.  Additionally,  this  work  demonstrated  the  ability  of  the  col¬ 
location  method  to  produce  a  more  simplified  formulation  than  the  Galerkin 
method.  The  successes  reported  in  this  paper  motivated  the  idea  of 

developing  a  numerical  method  based  on  collocation  for  solving  the  EPFBE. 

**  *  * 

To  insure  that  this  attempt  did  not  violate  the  second  criterion,  a 
special  trial  solution,  derived  using  the  concept  of  space-angle-synthesis, 
was  used. 

This  combination  of  trial  solution  and  numerical  method  allowed  the 
successful  approximation  of  the  EPFBE.  The  simplicity  of  the  trial 
solution  allowed  all  integrals  in  this  formulation  to  be  analytically 
evaluated  on  MACSYMA.  Attempts  at  solving  the  air-over-ground  problem 
demonstrated  a  numerical  problem  associated  with  the  definition  of  the 
odd-parity  fluence.  The  severity  of  this  problem  eliminated  the  use 
of  the  collocation  method  for  solving  the  EPFBE  when  a  standard  finite 
element  or  combination  finite  element  synthesis  trial  solution  was  used. 
Collocation  Method  and  Trial  Solution 

Collocation  is  classified  as  a  weighted  residual  method.  In  the 
previous  chapter  a  derivation  was  presented  to  illustrate  the  weighted 
residual  method.  From  this  derivation  the  following  equation  resulted. 

J  [  L  Z  fci)]ccxJcLt  *  O  (3.1) 

The  difference  between  the  collocation  method  and  Galerkin  method  is  in 
the  definition  of  the  weight  function, EC?;  .  As  previously  demonstrated 
this  weight  function  is  defined  in  the  following  manner  for  the  Galerkin 
method. 
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(3.2) 


M 


£(.xJ  -  IZkQ- 

i*  i 

For  the  collocation  method  the  weight  function  is  defined  as 
ECx)  ~  S(x-  x.) 


(3.3) 


where  is  the  Dirac  delta  function.  Substituting  Eq.  (3.3)  into 

Eq.  (3.1)  gives 

f  [L  Z  aL  Q.Cx)-  fcxJ ]<fc?-V.)c/x  -  O  (3.4) 


C*  1 


Since  the  Dirac  delta  function  is  defined  as 

ia-i.,  -  S  °  • 


(3.5) 


00  ,  *  =  X 


L 


x -  -x  )  cL*  - 


Eq.  (3.4)  becomes 


N1 

L  Z  a^QaxJ  s  T  cx„, 

i*  i 


(3.6) 


To  solve  for  the  a;  coefficients,  Eq.  (3.6)  must  be  specified  at  M 
points  in  the  domain  D.  Thus,  Eq.  (3.6)  can  be  written  as 

Z  Qi  L  Qi  csfj  =  -fcxj  (3-7) 


c=  i 


Eq.  (3.7)  is  a  general  expression  representing  the  collocation  method. 

For  the  EPFBE,  using  the  trial  solution  presented  in  Eq.  (2.16),  this 
expression  becomes 

M  N  \ 

Z  ^  a.1.  \  -  si-  V  V  SLlrJ  AjCAjih  (3.8) 

i ~l  Jr|  * 

-A  ^  -  '5gO\lA.jJ 
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The  coefficients  cuj  are  determined  by  solving  a  system  of  simultaneous 
equations  resulting  from  the  evaluation  of  Eq.  (3.8)  at  NXM  phase  space 
points. 

Comparing  Eq.  (3.8)  to  Eq.  (2.18)  demonstrates  the  simpl ication 
resulting  from  the  collocation  method.  Most  significant  is  the  absence 
of  the  weight  function  and  outer  integrals  that  caused  complications  in 
the  Galerkin  weak  form  equation.  Eliminating  these  two  complicating 
factors  allows  the  first  criterion  specified  earlier  to  be  met. 

The  collocation  method  does  introduce  some  new  complexities.  The 
most  obvious  is  that  the  vacuum  boundary  condition  is  no  longer  treated 
naturally.  This  condition  must  now  be  enforced  by  using 

Vcr, A)  =  -  sx-v  =  O,  A-rL  <  O  (3.9) 

as  presented  in  Chapter  I.  Unlike  the  Galerkin  method,  the  collocation 
method  does  not  generate  a  symmetric  global  matrix.  This  lack  of  sym¬ 
metry  and  positive  definiteness  increases  the  cost  of  solving  the  result¬ 
ing  matrix  and  eliminates  the  use  of  some  very  efficient  and  fast  linear 
equation  solving  algorithms.  Another  complication  results  from  the 
required  continuity  of  the  trial  solution.  In  the  Galerkin  weak  form 
equation,  the  functions  chosen  for  use  in  the  trial  solution  were  required 
to  be  continuous.  The  variational  nature  of  this  equation  eliminated 
continuity  constraints  on  the  derivatives  of  these  functions  and  thus 
allowed  the  use  of  piecewise  continuous  Lagrange  polynominals.  The 
collocation  method  requires  the  trial  solution  to  have  continuous  second 
derivatives  at  the  mesh  nodes.  The  simplification  gained  by  using  the 
collocation  method  forces  higher  continuity  requirements  on  the  trial 
solution. 
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A  trial  solution  requiring  continuity  in  the  second  derivative  may 
be  constructed  from  cubic  splines.  A  derivation  of  the  cubic  spline 
function  in  one  dimension  is  presented  in  Appendix  I  and  is  denoted  as 

.  A  tensor  product  consisting  of  cubic  splines  representing  both 
the  spatial  and  angular  variation  of  the  neutron  fluence  would  result 
in  a  global  matrix  of  large  bandwidth.  In  an  air-over-ground  problem 
both  the  spatial  and  angular  domain  demand  relatively  fine  zoning.  This 

zoning  requirement  increases  the  bandwidth  problem  in  the  global  matrix 

* 

and  increases  the  computational  cost.  To  reduce  the  bandwidth  a  synthesis 
function  was  used  instead  of  cubic  splines  to  represent  the  angular  varia¬ 
tion  of  the  neutron  fluence. 

An  angular  synthesis  function  was  devised  to  represent  the  angular 
variation  of  the  even-parity  fluence  and  meet  the  second  criterion 
described  previously.  This  function  can  be  represented  as 

Ai/JL.X)  =  a,  t  (  M/Wo  +  V'-ZZ/  cosCx-vJ2  (3-10) 

where  Uq  and  xq  are  the  coordinates  of  the  streaming  direction  with 
respect  to  the  point  source  in  the  particle  direction  coordinate  system. 

This  synthesis  function  was  devised  by  considering  the  expected  angular 
variation  of  the  neutron  fluence.  Appendix  G  describes  this  synthesis 
function. 

Eq.  (3.10)  has  properties  similar  to  the  ellipsoidal  synthesis  func¬ 
tion  used  in  Chapter  II.  Eq.  (3.10)  is  pitched  in  line  with  the  streaming 
direction  at  a  particular  spatial  location.  This  equation  can  also  be 
varied  to  represent  the  three  expected  angular  fluence  shapes  presented 
in  Fig.  6.  A  highly  anisotropic  flux  distribution  occurs  when  the 
coefficient  is  much  larger  than  the  a,  coefficient.  An  isotropic 
flux  distribution  occurs  when  ad  =  0. 
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The  primary  difference  between  the  synthesis  function  defined  in 
Eq.  (3.9)  and  the  ellipsoidal  synthesis  function  presented  in  Chapter 
II  is  the  meaning  of  the  coefficients.-  The  coefficients  in  the  ellipsoidal 
synthesis  function  expansion  allowed  the  blending  of  three  ellipsoidal 
shapes  to  determine  the  angular  approximation.  In  Eq.  (3.10)  the  coeffi¬ 
cients  directly  shape  the  synthesis  function  for  approximating  the  angular 
variation  of  the  even-parity  fluence. 

The  use  of  a  synthesis  function  to  represent  the  angular  dependence 
of  the  even-parity  fluence  reduces  the  bandwidth  of  the  resulting  global 
matrix.  This  reduction  occurs  since  the  synthesis  function  is  globally 
defined  over  the  entire  angular  domain,  thus  eliminating  the  need  for 
zoning.  One  globally  defined  function  replaces  all  the  tensor  product 
combinations  that  would  result  from  an  angular  trial  solution  defined 
over  a  discretized  domain.  Since  each  tensor  product  combination  occupies 
one  element  in  the  global  matrix,  the  globally  defined  synthesis  function 
can  significantly  reduce  the  bandwidth  of  this  matrix.  An  additional 
advantage  to  using  the  synthesis  function  defined  in  Eq.  (3.10)  is  that 
it  satisfies  the  reflective  boundary  condition.  The  final  form  of  the 


trial  solution  used  with  the  collocation  method  is 


'-I  M 

KfJLt,A)=  ^  (pj  8;CIJ  la,i:  *  Qii;  (’/UyU0*>T7-« 

t =  I  J r  1  v  J  V  J 


1 


(3.11) 


where  Qtpi  and  6JCH  are  cubic  splines  in  the  spatial  domain. 

The  combination  of  the  collocation  method  and  the  trial  solution 
presented  in  Eq.  (3.11)  satisfied  the  two  criteria  determined  from  Chapter 
II.  A  computer  program  was  written  that  used  this  combination  to  generate 
a  global  matrix  representation  of  the  EPFBE.  This  program  used  no  numer¬ 
ical  integration,  since  the  simplicity  of  the  trial  solution  and  the  col¬ 
location  method  allowed  all  integral  and  derivative  operations  to  be 
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performed  on  MACSYMA.  The  results  of  these  operations  were  translated 
into  FORTRAN  statements  by  MACSYMA,  thus  eliminating  the  possibility  of 
erroneously  transcribing  these  expressions. 

A  distributed  first-scatter  source  term  was  used  in  this  computer 
code.  The  angular  distribution  of  the  neutron  fluence  near  the  point 
source  is  extremely  anisotropic.  This  anisotropy  is  so  severe  that  the 
angular  synthesis  function  represented  by  Eq.  (3.10)  would  have  difficulty 

in  accurately  approximating  such  an  anisotropic  shape.  The  use  of  the 

.  - 

first-scatter  source  allows  the  scattered  and  unscattered  neutrons  at  a 
point  to  be  treated  separately.  Eq.  (1.6)  gives  at  r  the  angular  dis¬ 
tribution  of  first-scattered  neutrons  which  were  emitted  from  a  point 
source  located  at  »s  .  The  uncollided  neutrons  arriving  at  r  are  given 
by 

while  those  that  interact  either  through  absorption  or  scatter  are 
represented  by 


The  scattered  fluence  is  less  anisotropic  then  the  total  fluence  and  may 
be  represented  by  a  function  with  a  simple  angular  dependence.  In  the 
final  solution  the  uncoil ided  fluence  is  added  to  the  calculated  value 
of  the  neutron  fluence  in  the  streaming  direction. 

As  previously  mentioned,  the  collocation  method  requires  an  evalua¬ 
tion  of  Eq.  (3.8)  at  a  number  of  phase  space  points  equal  to  the  number 
of  coefficients  in  the  trial  solution.  The  trial  solution  presented  in 
Eq.  (3.11)  requires  Eq.  (3,8)  to  be  evaluated  in  two  distinct  angular 
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directions  at  each  spatial  location.  In  the  collocation  based  computer 
code,  these  angular  directions  were  always  selected  as  the  streaming 
direction  and  a  direction  perpendicular  to  the  streaming  direction.  The 
choice  of  these  two  directions  allowed  the  coefficients  to  be  defined 
in  the  most  meaningful  manner  based  on  their  role  in  properly  shaping  the 
angular  synthesis  function.  The  spatial  collocation  points  were  chosen 
as  the  nodes  of  the  mesh.  The  only  exception  occurred  along  the  mesh 
linep«o.  On  this  mesh  line  the  reflected  boundary  condition  must  be 
enforced.  Since  the  chosen  trial  solution  automatically  satisfied  this 
constraint,  spatial  collocation  points  were  chosen  at  a  value  off  slightly 
greater  than  zero.  Along  the  vacuum  boundaries  the  code  provided  the 
option  of  either  evaluating  Eq.  (3.9)  at  two  angular  directions  or  enforc¬ 
ing  the  EPFBE  in  one  angular  direction  and  the  vacuum  boundary  condition 
in  the  other  direction.  A  flow  diagram  of  the  collocation  computer  code 
appears  in  Appendix  H. 

The  structure  of  the  global  matrix  resulting  from  the  use  of  the 
selected  trial  solution  and  the  collocation  method  was  primarily  determined 
by  the  cubic  splines.  The  structure  of  this  matrix  for  the  spatial  mesh 
in  Fig.  7  is  illustrated  in  Fig.  8. 

Problems 

An  air-over-ground  problem  was  attempted  using  the  collocation 
method.  The  spatial  mesh  used  for  this  problem  was  identical  to  the  one 
illustrated  in  Fig.  7.  Collocation  points  were  selected  as  previously 
described  and  the  vacuum  boundary  condition  was  enforced  along  the  appro¬ 
priate  boundary  surfaces.  To  generate  a  solution  to  this  problem  the 
collocation  global  matrix  assembly  program  was  expanded.  This  expansion 
was  necessary  to  determine  the  coefficients  of  the  trial  solution. 
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Fig  8.  Global  Collocation  Matrix 
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analyze  the  resulting  global  matrix,  and  reconstruct  in  terms  of  the  trial 
solution  the  even-parity,  Boltzmann,  and  scalar  fluences.  Solution 
coefficients  were  determined  using  Gaussian  elimination  with  full  pivoting. 
The  input  to  this  linear  equation  solving  routine  was  the  collocation 
global  matrix  and  source  column  vector  computed  from  the  global  matrix 
assembly  program.  The  global  matrix  was  analyzed  by  determining  its  con¬ 
dition  number  which  was  calculated  as  the  ratio  of  the  largest  to  the 

smallest  eigenvalue.  (The  eigenvalues  were  determined  by  the  International 

*  * 

Mathematical  and  Statistical  Libraries  routine  entitled  LSVDF.)  Recon¬ 
struction  of  the  desired  fluences  was  accomplished  from  the  solution 
coefficients  and  basis  functions  of  the  trial  solution  for  the  spatial 
location  chosen.  The  reconstruction  program  used  many  of  the  same  sub¬ 
routines  that  were  developed  for  the  global  matrix  assembly  program.  The 
logic  of  this  system  of  programs  is  Illustrated  in  Fig.  9.  Each  subrou¬ 
tine  contained  within  the  global  matrix  assembly  program  and  fluence 
reconstruction  program  was  independently  checked  to  assure  that  its  logic, 
programming,  and  interface  with  the  driver  routines  was  correct. 

The  solution  generated  by  the  collocation  method  to  this  air-over- 
ground  problem  was  disappointing.  Most  significant  was  that  the  scalar 
fluence  was  negative  at  several  node  locations  in  the  spatial  mesh. 

Also  both  the  Boltzmann  and  even-parity  fluence  exhibited  negative  values 
for  several  angular  locations  at  each  spatial  node  point.  Negative  angular 
fluences  values  were  not  limited  to  spatial  node  points,  but  also  appeared 
at  off  node  spatial  locations  throughout  the  spatial  domain.  Those  scalar 
fluence  values  that  were  positive  did  not  make  physical  sense.  The  value 
of  the  scalar  fluence  was  often  larger  on  the  boundaries  then  near  the 
interior  of  the  spatial  mesh.  Also  the  condition  number  of  the  resulting 
global  matrix  was  10  . 


Fig  9.  Program  Log 


Several  means  were  attempted  to  improve  these  poor  results.  The 
meshing  of  the  spatial  domain  was  refined.  This  refinement  increased 
the  number  of  nodes  from  sixteen  to  thirty-six,  thus  increasing  the  size 
of  the  global  matrix  to  a  72  x  72  system.  This  refinement  produced  no 
change  in  the  previously  observed  behavior  of  the  fluences  and  increased 

O 

the  condition  number  to  10  . 

The  method  of  selecting  collocation  points  was  altered.  Spatial  col¬ 
location  points  were  selected  at  locations  other  than  at  the  spatial  nodes 

♦ 

and  angular  collocation  points  were  selected  away  from  the  streaming 
direction.  Again  the  behavior  of  the  fluence  remained  the  same  and  the 
conditon  number  did  not  change  appreciably  from  the  original  value  of  105. 
In  addition  a  new  difficulty  was  encountered.  This  difficulty  involved 
selecting  spatial  collocation  points  that  were  not  located  on  either  the 
vertical  or  horizontal  mesh  lines.  Distributing  these  points  evenly 
around  the  spatial  mesh  proved  impossible  and  forced  a  random  type  of 
selection.  The  optimum  way  of  picking  these  points  from  the  unlimited 
number  of  combinations  available  was  not  clear.  This  random  selection 
process  was  abandoned,  since  the  resulting  solution  exhibited  the  same 
problems  originally  encountered. 

Another  attempt  at  resolving  these  solution  problems  used  the  option 
mentioned  earlier  of  collocating  either  the  EPFBE  or  vacuum  boundary 
condition  on  the  appropriate  boundary  surface.  Two  different  approaches 
were  tried.  The  first  approach  allowed  the  vacuum  boundary  condition  to 
be  used  for  both  angular  collocation  points  at  a  boundary  spatial  node. 

The  second  approach  divided  the  two  angular  collocation  evaluations 
between  the  EPFBE  and  vacuum  boundary  condition.  In  both  cases  the 
behavior  of  the  scalar  fluence  and  value  of  the  condition  number  remained 


51 


effectively  the  same  as  the  original  result.  At  this  point  it  was  becom¬ 
ing  increasingly  clear  that  a  fundamental  problem  existed  with  the  colloca¬ 
tion  method.  To  uncover  this  problem  a  simple  diagnostic  problem  was 
used. 

This  diagnostic  problem  assumed  a  uniform  isotropic  neutron  source, 
isotropic  scatter,  and  isotropic  boundary  conditions.  Its  solution,  as 
verified  by  a  DOT  3.5  calculation,  is  a  flat  spatial  distribution  isotro¬ 
pically  distributed  in  angle.  The  simplicity  of  this  problem  allowed 
several  different  formulations  of  the  EPFBE  to  be  evaluated.  These'*various 
formulations  differed  in  complexity  and  all  had  the  capacity  of  generating 
an  accurate  solution  to  this  problem. 

The  simplest  formulation  was  derived  by  using  the  isotropy  of  the 
known  solution.  This  isotropy  allowed  all  terms  in  the  EPFBE  related  to 
the  odd-parity  fluence  to  be  eliminated.  The  resulting  formulation  was 

Gj  V(t.A)*  (3-14) 

The  collocation  global  matrix  assembly  program  was  modified  to  solve  Eq. 
(3.14).  This  modification  included  replacing  the  first-scatter  source 
term  and  vacuum  boundary  condition  with  those  associated  with  the  diagnos¬ 
tic  problem.  The  source  term  and  boundary  condition  for  the  diagnostic 
problem  were  respectively  defined  as 

'  Wn  (  <rT(?)  (3.15) 

T =  hT?  for  F*  on  the  boundary,  (3.16) 

The  collocation  method  generated  the  correct  solution  to  the  diagnostic 
problem  using  this  formulation  of  the  EPFBE.  This  result  was  encouraging. 
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since  it  confirmed  the  logic  and  programming  accuracy  of  the  subrountines 
in  the  global  matrix  assembly  program  and  reconstruction  program. 

A  more  complicated  formulation  was  attempted  next.  This  formulation 
was 

“  A  V  kp.  A  V  Vef'A)  *  Gj  -  SyCfrA)  (3.17) 


where  S*<*,*>  and  the  boundary  conditions  were  the  same  as  the  previous 
problem.  For  this  formulation  the  collocation  method  did  not  generate  a 
flat  scalar  fluence  distribution.  The  boundary  scalar  fluences  were  cor¬ 
rect,  but  the  interior  mesh  values  varied  by  as  much  as  20  percent.  For 
this  formulation  all  the  scalar  fluences  were  positive.  The  same  formula¬ 
tion  was  solved  again  using  a  more  refined  spatial  mesh.  The  solution 
was  the  same,  indicating  the  original  mesh  was  adequate. 

The  formulation  presented  in  Eq.  (3.14)  was  next  used  with  the  fol¬ 
lowing  form  of  the  diagnostic  problem  boundary  condition 


YdAl  *  Kcr,*)- 


HTt 


(3.18) 


Eq.  (3.18)  represents  a  statement  of  the  diagnostic  problem  boundary 
condition  in  terms  of  the  Boltzmann  fluence.  The  collocation  method  must 
be  able  to  handle  this  type  of  boundary  condition  formulation,  since  in 
all  transport  problems  the  physical  constraints  on  the  boundaries  are 
specified  in  terms  of  this  fluence.  In  the  collocation  computer  program 
XrtAi  was  replaced  by  the  odd-parity  fluence  transformation  Eq.  (A. 20). 

The  solution  generated  by  the  collocation  method  for  this  problem  was  again 
incorrect.  The  interior  scalar  fluences  were  correct;  however,  the 
boundary  scalar  fluences  varied  by  120  percent. 

The  final  formulation  used  Eq.  (3.17)  with  the  boundary  condition  pre¬ 
sented  in  Eq.  (3.18).  The  collocation  method  generated  a  solution  for 
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this  formulation  that  was  utter  nonsense.  This  solution  had  several  nega¬ 
tive  fluence  values  both  in  the  interior  of  the  mesh  and  on  the  boundary. 

The  scalar  fluences  that  were  positive  were  not  equal  and  could  not  have 
generated  a  flat  spatial  distribution.  ' 

Problem  Analysis 

The  poor  results  generated  by  the  collocation  method  can  best  be 

explained  by  using  a  1-D  plane  geometry  formulation  of  the  diagnostic 

problem.  In  this  simplified  formulation  a  suitable  synthesis  trial  solu- 

*  * 

tion  for  the  standard  Boltzmann  equation  is 

@LCx)(a,i  tCL2iM)  (3-19) 

where  @L(x)is  a  cubic  spline  basis  function.  Using  the  defintion  of  the 
even-parity  fluence  with  Eq.  (3.19)  generates  the  following  trial  solution 

Y 2  a,-£tcx)  (3*20) 

c—  I 

Applying  the  odd-parity  fluence  transformation  presented  in  Eq.  (A. 20)  to 
the  trial  solution  of  Eq.  (3.20)  results  in 

X(x,u)  s  ~  Z  a,i  £oT>  &  (%(*>)  (3-21> 

L-i  ' 


for  a  PQ  scatter  process.  The  boundary  condition  of  the  diagnostic 
problem  in  terms  of  'Pcx.kI  and  can  be  written  as 


T* 


-L.Cx.At)  =  YcK'V./  *  Xcx.UL)  -  Hn 


(3.22) 


Using  Eqs.  (3.20)  and  (3.21),  Eq.  (3.22)  can  be  rewritten  as 
CX,IL)  ~  a,i  (  Q(Xj  ~  OfU)  ~ctx  (Q(x>)) 


(3.23) 
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By  the  definition  of  the  collocation  method  Eq.  (3.23)  can  be  satisfied 
for  one  selected  angular  collocation  point  at  a  particular  spatial  node. 

Eq.  (3.23)  can  not  be  properly  constrained  to  generate  the  isotropic 
boundary  condition.  In  any  other  angular  direction  other  than  the  selected 

I 

collocation  direction  the  value  of  the  boundary  condition  is  not  sn  . 

The  odd-parity  fluence  transformation  has  forced  the  Boltzmann  fluence  to 
have  an  anisotropic  component.  This  can  be  easily  seen  by  the  presence  of 


the/*  variable  in  Eq.  (3.23).  Ideally  for  the  diagnostic  boundary  condition 


the  following  relation  should  hold: 


~Xcx,n)=  Za.i  O 

is  i  * 


(3.24) 


however;  this  would  require  either 

(1)  all  the  a.j'jto  be  zero  in  which  case  %<,*</*  O  at  all  phase 
space  points. 

(2)  A  (t c  O  ,  which  is  only  true  for  the  spline  centered  at 

a  node  and  since  splines  overlap  in  mesh  intervals  this  is  not 
in  general  true  (  Xe-  boundary  node  location). 

(3)  cla  &(&**}[+ JJ+* &  &  (&*M)ljp 


which  is  not  enforced  in  the  collocation  set  of  equations. 

The  conclusion  is  that  the  selected  synthesis  trial  solution  is  not 
transformable  through  the  odd- parity  fluence  transformation  into  a  trial 
solution  that  accurately  represents  the  anisotropy  of  the  fluence. 

The  odd-parity  fluence  transformation  is  also  present  in  the  even- 
parity  equation.  As  demonstrated  in  Appendix  A  the  odd-parity  fluence  is 
contained  in  the  first  term  of  this  equation.  If  the  even-parity  equation 
is  being  derived  for  the  diagnostic  problem,  this  first  term  would  appear 
as 
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M. 


£L- 

CLX 


(  X(*,iu) 


(3.25) 


To  eliminate  the  X  variable  from  the  even-parity  equation  the  odd-parity 
fluence  transformation  is  used.  This  is  accomplished  by  substituting  the 
value  of  5CC from  Eq.  (3.24)  into  Eq.  (3.25),  thus  generating 


A 


(-  A _ 

Uk  (  ctcki 


(3.26) 


Simplifying  Eq. 


(3.27) 


The  full  even-parity  equation  for  the  1-D  plane  geometry  diagnostic  prob¬ 
lem  is 

mJ~  0J4<x)f VtX'U'j  c(p.‘  ~  (3.28) 

A* 

Substituting  the  trial  solution  of  Eq.  (3.20)  into  Eq.  (3.28)  yields 

zE  a.;  (  oyu)  c[$t~  *  (o-ji*)  -  o^c*))  &cc*>')z  "477  (3.29) 

i 


X.JULI 


crT(x) 


CUT 


cr.CxI 


fcx, 


Obviously  to  generate  a  flat  isotropic  Boltzmann  fluence,  the  first  term 

of  Eq.  (3.28)  must  be  zero,  since  its  origin  is  the  difference  between 

./ul  el*  Pit*) 

the  leakage  in  the/*  and-/*  angular  directions.  If  <rT(*>  ~~o7r'  could 
be  made  equal  to  zero  then  Eq.  (3.29)  would  be 


w  / 

Z  ( <rTu)  -  <%*(.*))  &c*\  =  (3.30) 

IM  c 


Z  a,i 


= 


) _ _ 

4n  (cf(t)-fy  (.*)) 


(3.31) 


which  generates  the  correct  solution  to  the  diagnostic  problem.  In  the 
collocation  method  there  is  no  means  to  set  this  term  equal  to  zero. 


56 


Just  as  in  the  boundary  condition,  the  odd-parity  fluence  transformation 
has  forced  an  incorrect  anisotropy  into  the  solution. 

To  determine  tne  form  of  the  anisotropy  resulting  from  the  odd-parity 
fluence  transformation,  an  analysis  was  performed.  This  analysis  was 
based  on  an  interpolation  procedure  similar  to  the  one  used  to  apply  the 
Boltzmann  fluence  boundary  condition  to  the  even-parity  equation.  In 
this  procedure  three  different  trial  solutions  were  used  to  represent 
the  even-parity  fluence.  These  trial  solutions  were  selected  because 
they  illustrated  specific  problems  related  to  the  odd-parity  fluence  trans¬ 
formation.  For  each  trial  solution  a  simple  spatial  mesh  was  constructed 
and  where  necessary  an  angular  mesh.  All  results  were  generated  using 
the  same  sequence  of  programs  that  were  used  in  solving  the  four  versions 
of  the  diagnostic  problem.  Modifications  to  these  programs  were  made  to 
account  for  the  different  interpolating  functions.  All  angular  plots 
presented  are  typical  of  those  seen  at  all  the  nodes  in  the  spatial  mesh. 

The  first  problem  used  the  trial  solution  presented  in  Eq .  (3.19)  and 
solved  the  interpolation  problem  of  Eq.  (3.22)  with  the  1/4  rr  repl aced  by 
100.  Fig.  10  illustrates  the  results  of  this  problem  for  various  values 
of  the  total  macroscopic  cross  section  (ojz*;) .  This  interpolation  problem 
demonstrates  graphically  the  explanation  that  was  presented  using  Eqs. 
(3.19)  thru  (3.31).  The  spatial  mesh  for  this  problem  had  nodes  at  0, 

10,  15,  and  20  cm.  The  angular  collocation  point  at  each  of  these  spatial 
nodes  was  taken  at  u  =  1  which  corresponds  to  the  streaming  direction  in 
a  1-D  plane  geometry  problem. 

-4 

In  Fig.  10  <rT  was  varied  from  infinity  to  10  /cm  thus  varying 


the  effect  of  the  odd-parity  fluence  transformation.  For«rr=  <*>  an 
isotropic  fluence  distribution  across  the  angular  domain  was  generated. 


r 


This  result  was  analogous  to  the  first  version  of  the  diagnostic  problem 
that  generated  acceptable  results.  As  illustrated  by  the  other  curves  in 
Fig.  10,  reducing  the  value  of  crr  increases  the  weight  of  the  odd-parity 
fluence  transformation  in  determining  the  angular  distribution  of  the 
Boltzmann  fluence.  The  effect  of  this  transformation  is  most  clearly 
demonstrated  by  the  odd  functional  nature  of  the  resulting  Boltzmann 
angular  fluence  distribution  for  crr  =  \o'H  /cm. 

For  the  trial  solution  presented  in  Eq.  (3.20),  the  resulting  aniso¬ 
tropy  is  fixed.  The  collocation  method  can  alter  the  magnitude  of  the 
anisotropy,  but  can  not  change  its  relative  angular  distribution.  Thus 
the  anisotropy  demonstrated  in  Fig.  10  will  exist  at  all  spatial  locations. 
This  is  obviously  unsatisfactory  since  the  anisotropy  in  most  transport 
problems  (especially  in  an  air-over-ground  problem)  varies  from  one 
spatial  location  to  another. 

In  the  second  interpolation  problem  a  pure  spline  basis  is  used  as 
a  trial  solution  for  the  even-parity  fluence.  This  trial  solution  can 
be  represented  mathematically  as 


l,j?  a‘J 


(3.32) 


where  £?<yjand  ^Y^/are  cubic  splines  defined  respectively  over  the 
spatial  and  angular  domain.  From  this  trial  solution  the  Boltzmann  flu¬ 
ence  is  formed  for  a  1-D  plane  geometry  and  appears  as 


-  2]  7z,  o.ij  (friij&iAii  ~  W  ^  )  (3 


33) 


Eq.  (3.33)  was  interpolated  across  a  spatial  mesh  with  three  nodes  located 
at  1.,  10.,  and  20.  cm  and  an  angular  mesh  with  three  nodes  located  at 
-1,  0,  and  1.  The  spatial  and  angular  nodes  also  served  as  collocation 
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coordinates,  where  the  Boltzmann  angular  fluence  (Eq.  (3.33))  was  forced 
to  equal  100.  Again  with  <rT  set  equal  to  infinity  a  flat  isotropic 
angular  distribution  was  generated  across  the  spatial  domain.  In  Fig. 

11  and  Fig.  12  the  value  of  o-T  i s  continuously  decreased.  The  curves  in 
these  figures  demonstrate  the  same  trend  seen  previously  with  the  1-D 
synthesis  trial  solution.  As  <rT  decreases  in  value  the  dominance  of 
the  odd-parity  fluence  transformation  increases. 

Unlike  the  synthesis  trial  solution  the  pure  spline  trial  solution 
can  alter  the  anisotropy.  This  alteration  can  be  achieved  by  more  finely 
discretizing  the  angular  domain  and  thus  forcing  the  solution  to  be  con¬ 
strained  at  more  angular  nodes.  This  approach  generates  two  disadvantages 
First,  the  behavior  of  the  trial  solution  in  Eq.  (3.33)  would  still  be 
dominated  by  the  odd-parity  fluence  transformation  between  the  nodes  for 
small  values  of  erT  .  This  dominance  forces  an  artifical  anisotropy 
between  the  nodes  that  in  an  exact  solution  would  not  exist.  This  arti¬ 
fical  anisotropy  would  affect  the  evaluation  of  the  boundary  condition 
and  first  term  of  the  even-parity  equation  in  a  manner  similar  to  that 
observed  in  the  diagnostic  problem.  The  degree  to  which  this  effect 
would  alter  the  final  solution  is  directly  related  to  the  number  of  nodes 
used  in  the  angular  domain.  Increasing  the  number  of  nodes  leads  to 
the  second  disadvantage.  This  disadvantage  addresses  the  size  of  the 
resulting  global  matrix  that  would  be  generated  using  a  pure  spline 
basis.  The  total  number  of  coefficients  that  must  be  solved  for  to  obtain 
a  solution  to  the  EPFBE  is  equal  to  the  number  of  unknown  coefficients  in 
the  trial  solution.  The  trial  solution  presented  in  Eq.  (3.33)  would 
result  in  a  global  matrix  that  has  a  row  and  column  dimension  equal  to 
NXM.  Increasing  the  number  of  nodes  by  A  in  the  angular  domain  results 
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Fig  11.  Interpolation  Problem  2  (1  of  2) 


in  an  increase  in  the  dimensionality  of  the  global  matrix  of  M  .  For 
an  air-over-ground  problem  a  suitable  pure  spline  trial  solution. for  the 
even-parity  fluence  is 


Z  al/L  6  (fi  Sjtii  (3.34) 

J 


Refining  the  angular  mesh  a  factor  of  A  in  the /x  and  X  coordinates  would 
result  in  a  (otpjt  increase  in  the  dimensionality  of  the  glob¬ 

al  matrix.  The  acceptable  value  of  A  needed  to  offset  the  artific^l 
anisotropy  effect  of  the  odd-parity  fluence  transformation  is  not  known, 
but  based  on  the  results  of  the  diagnostic  problem  and  presented  figures 
would  appear  to  be  quite  large.  A  large  A  would  result  in  a  large 
global  matrix.  This  matrix  would  be  nonsymmetric  and  costly  to  compute 
and  solve,  thus  eliminating  the  pure  spline  basis  as  an  efficient  means 
to  solve  the  even-parity  equation,  especially  for  an  air-over-ground 
problem. 

The  third  trial  solution  used  to  demonstrate  the  anisotropy  effect 
of  the  odd-parity  fluence  transformation  was  the  one  presented  in  Eq.  (3.11). 
For  this  trial  solution  an  interpolation  problem  was  not  formulated.  A 
more  simple  and  direct  means  of  demonstrating  the  effect  of  the  odd-parity 
fluence  transformation  was  found.  A  two-dimensional  spatial  mesh  in 
the  coordinates  p  and  i  was  constructed.  The  nodes  in  both  the  p  and 
l  direction  were  located  at  1.,  10.,  and  20.  cm.  At  the  coordinate  loca¬ 
tion  (10.,  10.)  the  Boltzmann  angular  fluence  distribution  was  calculated. 
This  calculation  was  performed  by  the  fluence  reconstruction  program. 

Input  to  this  program  included  the  spatial  mesh,  solution  coefficients, 
c\(.?)  ,  and  the  location  of  the  point  source  that  determines  the  pitch 
of  the  trial  solution  at  the  various  spatial  locations.  All  solution 
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mum 


coefficients  were  set  equal  to  one.  The  point  source  was  located  at  the 
spatial  coordinates  (0.,  10.),  thus  making  the  pitch  axis  of  the  trial 
solution  parallel  with  the  p  axis. 

Figs.  13  and  14  illustrate  the  results  of  the  fluence  reconstruction 
program.  Figs.  13a  and  13b  demonstrate  the  angular  distribution  of  the 
even-parity  fluence  for  x  =  0  and  77  respectively.  Figs.  13c  and  13d 
show  the  form  of  the  odd-parity  fluence  derived  from  the  trial  solution 
of  Eq.  (3.11)  for  various  x  locations.  Figs.  14a  thru  14d  are  plots  of 
the  Boltzmann  fluence  for  various  x  locations  and  values  of  <rr  .  Tlie 
significant  thing  to  note  from  these  plots  is  the  anisotropy  that  is 
forced  on  the  Boltzmann  fluence  by  the  odd-parity  fluence  transformation. 
A  polar  plot  of  Figs.  14a  and  14b  is  presented  in  Fig.  15.  This  polar 
plot  represents  the  angular  synthesis  function  for  the  Boltzmann  fluence 
as  derived  from  the  even-parity  fluence  angular  synthesis  function  and 
odd-parity  fluence  transformation.  This  function  is  not  representative 
of  the  angular  distribution  expected  from  either  the  diagnostic  problem 
or  an  air-over-ground  problem.  This  poor  representation  is  caused  by 
the  inadequacy  of  the  selected  even-parity  angular  synthesis  function  to 
be  transformed  into  an  accurate  representation  of  the  anisotropic  com¬ 
ponent  of  the  Boltzmann  fluence.  As  evidenced  by  Figs.  14c  and  14d  the 
derived  angular  synthesis  function  becomes  worse  as  the  dominance  of 
the  odd-parity  transformation  increases.  The  negative  values  seen  in 
Fig.  14c  could  have  contributed  to  the  negative  fluence  values  calculated 
in  the  diagnostic  and  air-over-ground  problem. 

The  conclusion  of  this  chapter  is  that  the  collocation  method  is 
not  suitable  for  solving  the  EPFBE.  The  main  failure  of  this  technique 
lies  in  its  point-wise  definition.  This  definition  forces  the  selected 
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c.  cr-i  x--o  ar-i  x=  V 


Fig  13.  Even  and  Odd-Parity  Fluence 
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trial  solution  to  accurately  represent  the  even-parity  fluence  and  be 
transformable  into  an  accurate  representation  of  the  odd-parity  fluence. 
As  seen  in  this  chapter  the  standard  finite  element  trial  solution  (pure 
spline  basis)  and  combination  finite  element  synthesis  function  trial 
solution  do  not  meet  this  criterion.  The  Galerkin  method  reduces  these 
stringent  requirements  on  the  trial  solution  by  allowing  the  even-parity 
equation  to  hold  in  an  integral  sense  rather  than  a  point-wise  sense. 

The  integration  operations  associated  with  the  Galerkin  weak  form  equa¬ 
tion  act  as  a  smoothing  apparatus  and  shapes  the  trial  solution  to 'meet 
the  phase  space  volume  requirements  imposed  by  the  weak  form  equation. 
Ironically  the  very  thing  that  caused  the  Galerkin  method  to  fail,  the 
nested  integrals,  is  also  the  reason  a  finite  element  method  can  success¬ 
fully  generate  solutions  to  the  EPFBE. 
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IV.  Multigroup  Considerations 


Background 

The  energy  dependence  of  neutral  particle  fluences  must  be  considered 
in  an  air-over-ground  problem.  The  multigroup  method  (Ref.  1)  is  used  to 
approximate  this  dependence.  In  multi  group  theory  the  definition  of 
various  cross  sections  are  modified  and  new  cross  sections  are  defined. 
Multigroup  cross  sections  have  been  previously  presented  in  relation  to 
Eq.  (1.3)  for  the  regular  Boltzmann  equation.  Similar  definitions  result 
from  the  EPFBE  when  formulated  in  a  multigroup  structure.  The  EPFBE 
written  for  group  h  is 


-  A-  A  vYct?)  +  G*  =  <5r(vt)-Av  Kk  S*£*,*J 


where 

t  (/; A)  =  the  even-parity  fluence  for  group  h 
and  Gh  are  equivalent  to  the  Ku  and  Gg  operators  written  for  group  h 
and  are  defined  as 


k;  fa,  * 

Gl fa, * 


<j*cn  (  f  (?) 

iXAiitA'idx ) 

(4.2) 

af(rj  f(A) 

f  0-^  (A-AJ  f  (A')clA‘ 

(4.3) 

where 

oft?)  =  macroscopic  total  cross  section  for  group  h 

J  =  in  group  even-parity  differential  scattering  cross  section 


O  — .  A 

(  oft?) 


CO 

-  I 

Ji-odd. 


Cl\?' 

af(fi)  -  ofc?) 
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=  j£th  order  Legendre  coefficient  for  group  h) 


(  P^(A'A) 


=  jtth  order  Legendre  polynomial) 
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or(^Aj  and  or  £>■,*.  /  are  equivalent  to  the  even  and  odd-parity 
source  terms  written. for  group  h.  These  terms  differ  from  the  monoenergetic 
forms  of  these  source  terms  in  that  the  down-scatter  from  higher  energy 
groups  must  be  included  as  dictated  by  particle  conservation.  The  even 


and  odd-parity  multigroup  source  terms  can  be  defined  respectively  as 


^Srfr.AK  Se\r,nJ  +  f  Z  Oltu  t,  a  a  J  Y  ( Met* 

1  A/JJ  ij=l  ^ 

r  K.i 

ST  (?,  AJ  “  (r,sl)  J  ^  Cf  ,y\.A.  )"jC  (<\)gIa- 


(4.4) 

(4.5) 


where 

=  group  h  even-parity  source  term 
o0  ii-./y.]  =  group  h  odd-pan ty  source  term 
(f,A'A  )  =  even-parity  down  scatter  cross  section  from  group  g  to 
group  h 

=  odd-parity  down-scatter  cross  section  from  group  g  to 
group  h 

I  <-r  A.)  =  even-parity  fluence  in  group  g 

A  =  odd-parity  fluence  in  group  g 

Multigroup  Analysis 

A  multigroup  solution  of  Eq.  (4.1)  was  accomplished  by  Kaper,  Leaf, 
and  Lindeman  (Ref.  27)  using  the  variational  formulation  of  this  equation. 
These  authors  assumed  that  the  scattering  was  isotropic.  In  an  air-over¬ 
ground  problem  this  assumption  cannot  be  made,  since  scattering  in  air 
is  extremely  anisotropic.  When  the  scattering  must  be  modeled  as  aniso¬ 
tropic,  the  left  hand  side  of  the  weak  form  of  the  EPFBE  becomes  extremely 
complex  as  evidenced  in  Chapter  II.  When  the  multigroup  method  is  applied 
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to  either  the  EPFBE  or  its  weak  form  the  source  terms  increase  substan¬ 
tially  in  complexity.  The  multigroup  method  does  not  significantly  affect 
the  left-hand  side  of  either  formulation. 

To  illustrate  this  increased  complexity  a  comparison  is  made  between 
using  an  isotropic  and  anisotropic  scattering  representation  in  the  source 
terms  of  the  multi  group  EPFBE.  For  an  isotropic  scattering  process  Eqs. 
(4.4)  and  (4.5)  can  be  rewritten  as 


=  SeCr,A)*  Z  OlV  f  V  c//i' 

*  Jut j 

Sr(r,A.)  =  S0C?,Si) 


(4.6) 


(4.7) 


The  integral  term  !n  Eq.  (4.5)  is  only  defined  for  odd  Legendre  coefficients 
and  thus  for  an  isotropic  scattering  process  equals  zero.  Eq.  (4.6)  is  analo¬ 
gous  to  Eqs.  (2.28)  and  (2.29)  in  terms  of  simplicity  and  once  again  allows 
the  integral  in  Eq.  (4.6)  to  be  evaluated  once  to  any  desired  accuracy 
and  stored  for  future  use.  The  simplifying  assumption  of  isotropic  scat¬ 
ter  played  an  essential  role  in  the  successful  results  generated  by  the 
previously  mentioned  authors. 

Anisotropic  scattering  forces  the  use  of  the  odd-parity  fluence  term 
in  a  multigroup  formulation.  This  term  appears  in  the  odd-parity  source 
term  as  illustrated  in  Eq.  (4.5).  Substituting  the  definition  of  the  odd- 
parity  fluence  into  Eq.  (4.5)  yields 

■5r  Cr)Aj  ~  So  i.rA)  +  j'  Z  ^ 


The  complexity  of  Eq.  (4.8)  can  best  be  illustrated  by  assuming  a  multi¬ 
group  problem  with  a  groups.  For  the  first  group  the  odd-parity  source 
terni  can  be  written  as 
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The  second  group  odd-parity  source  term  must  include  a  down-scatter  term 
to  account  for  the  neutral  particles  that  scattered  from  group  1  to 


group  2. 


•X  C^a  I  ~  Secf,<xl  *  j  (Ts°n  (fA-A.)  JCcr,*)  cl  a.' 

*  /Liw 


(4.10) 


but  as  previously  described 


X  =  K,  (  5r(?,A')  "  A- 7  t(r,A )  ) 


(4.11) 


and  thus  Eq.  (4.10)  can  be  written  as 


f  cs'tfX/ijfk,  v  Vo*,*))} clA’  (4‘12) 

An 

The  complexity  of  the  multigroup  method  becomes  evident  in  formulating 
the  odd-parity  source  term  for  group  3.  This  source  term  can  be  written 
as 

'  ■'HU 


Expanding  the  second  term  in  the  integral  results  in 

/  ^(^a'a;  [  Kt  (  Stc?a')  -  a  v  fa-))}  dr 


(4.14) 


Substituting  Eq.  (4.12)  into  the  first  term  of  Eq.  (4.14)  generates 
/  o£a Crt A  a  )  [ S}»U-j +[&(**■*  )[K*(50L?a~J -A"  vf( ?,Ar})}ctx]  (4.15) 


If  all  the  terms  resulting  from  the  Ku  operator  were  written  out  in  full 
form,  one  term  in  Eq.  (4.15)  would  appear  as 


«cW a£tit,£&) j ) rcr'i?, AT"  A.-) SjAr‘id*rtixr‘c/A',ctA • (4 ‘16) 

An  An  An  -An 

In  a  similar  manner  the  term  in  Eq.  (4. .15)  containing  -*’■ v  ^ also 
becomes  embedded  in  a  set  of  nested  integrals.  When  the  operator  is 
fully  expanded  one  term  of  this  expansion  would  appear  as 

To*. ^ %5*-i (4-17) 

J"  A"  An  M-circU- 

The  nesting  problem  for  both  Eqs.  (4.16)  and  (4.17)  becomes  worse  when 

m  * 

the  number  of  energy  groups  increases.  The  problems  involved  with  eval¬ 
uating  the  integrals  in  Eq.  (4.17)  are  similar  to  those  mentioned  in 
Chapter  II  for  both  a  pure  spline  and  spline  synthesis  trial  solution. 

The  complexity  illustrated  in  Eqs.  (4.16)  and  (4.17)  increases  when 
the  odd-parity  source  is  placed  in  the  Galerkin  or  collocation  formulations. 
For  the  Galerkin  formulation  the  odd-parity  source  term  appears  in  the 


following  term 


L/  <  A-  vccf  A),  /(*  >  otAoCf 

-  V  •  r 


(4.18) 


In  the  collocation  method  it  appears  as 
/V  V  fC  XCrA.) 


(4.19) 


The  complexity  of  doing  multigroup  calculations  using  the  EPFBE  does 
not  exist  for  the  first-order  Boltzmann  equation.  The  essential  difference 
is  that  the  first-order  Boltzmann  "luence  is  not  defined  in  terms  of  a 
source.  The  odd-parity  f luence  is  dependent  on  the  odd-parity  source  term 
due  to  the  requirement  of  defining  this  f luence  in  terms  of  the  even- 
parity  fluence  during  the  derivation  of  the  EPFBE.  The  complexity  intro¬ 
duced  by  applying  the  multigroup  method  to  the  EPFBE  is  analogous  to  the 
complexity  that  evolved  in  the  Galerkin  method. 
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This  chapter  has  demonstrated  the  problems  involved  with  using  the 
EPFBE  in  a  mul tigroup  scheme.  These  problems  must  be  solved  before  this 
equation  can  be  used  to  generate  accurate  and  cost  effective  solutions 
to  an  air-over-ground  problem. 


V.  Summary,  Conclusions  and  Recommendations 
Summary  and  Conclusions 

The  purpose  of  this  research  was.  to  determine  the  feasibility  of 
using  the  even-parity  equation  to  calculate  the  neutral  particle  fluence 
distribution  for  an  air-over-ground  problem.  In  pursuing  this  research 
several  different  approaches  were  attempted.  Each  of  these  approaches 
involved  an  application  of  a  standard  finite  element  technique  to  an 
appropriate  form  of  the  EPFBE.  The  only  deviation  from  this  stand^cd 
approach  involved  the  use  of  the  synthesis  method.  This  deviation  can 
not  be  considered  radical  since  numerous  authors  reported  its  successful 
use  in  solving  transport  problems.  The  choice  of  the  EPFBE  was  motivated 
by  its  mathematical  properties  that  potentially  eliminated  ray  effects. 
Combining  the  finite  element  and  synthesis  method  to  solve  the  EPFBE 
appeared  to  be  reasonable. 

Two  basic  problems  related  to  the  EPFBE  evolved.  The  first  problem 
involved  the  complexity  that  occurred  when  the  weak  form  of  the  EPFBE  was 
applied  to  the  air-over-ground  problem.  The  two-dimensional  cylindrical 
geometry  required  by  this  problem  coupled  with  the  need  for  a  P3  scatter¬ 
ing  representation  contributed  heavily  to  this  complexity.  Another 
contributor  was  the  finite  element  method.  The  discretized  angular  and 
spatial  domain  used  with  this  method  and  the  linear  Lagrange  polynomial 
trial  solution  defined  over  these  domains  eliminated  the  use  of  analytical 
methods  for  evaluating  the  weak  form  equation  integrals.  Numerical  inte¬ 
gration  of  this  trial  solution  over  the  same  domains  was  extremely 
inefficient  and  would  have  required  a  high  order  quadrature  set  to 
generate  acceptable  results.  The  cost  of  numerically  evaluating  the  weak 
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form  integrals  to  the  necessary  accuracy  proved  prohibitive  for  a  full 
scale  air-over-ground  problem. 

This  integral  problem  originated  from  the  choice  of  the  Galerkin 
method.  Angular  integrals  were  already  present  in  the  EPFBE  because  of 
the  Ku  and  Gg  operators.  The  definition  of  the  Galerkin  method  added 
additional  angular  and  spatial  integrals  that  resulted  in  the  nesting  of 
up  to  six  integrals  in  some  terms  of  the  weak  form  equation.  Other  finite 
element  methods  such  as  least  squares  would  have  generated  the  same^ 
nesting  problems  and  complex  integrands  as  the  Galerkin  method. 

The  synthesis  method  was  not  successful  with  the  weak  form  of  the 
EPFBE  due  to  the  choice  of  the  function  representing  the  angular  distribu¬ 
tion.  This  synthesis  function  was  not  analytically  integrable  and  was 
complex  enough  to  require  the  use  of  higher  order  quadrature  sets. 
Additionally,  the  form  of  the  function  did  not  allow  the  pitch  parameters 
of  the  ellipsoid,  which  are  spatially  dependent,  to  be  separated  from 
the  angularly  dependent  variables.  This  fault  eliminated  the  possibility 
of  separating  the  angularly  and  spatially  dependent  terms  in  the  weak 
form  equation  and  forced  all  integral  evaluations  to  be  coupled.  Such  a 
coupling  required  all  integrals  to  be  reevaluated  whenever  the  spatial 
mesh  was  altered. 

A  conclusion  of  this  research  is  that  the  use  of  a  pure  linear 
Lagrange  trial  solution  coupled  with  any  of  the  integrally  defined  finite 
element  solution  techniques  does  not  provide  a  viable  means  for  solving 
the  EPFBE  formulated  for  an  air-over-ground  problem.  This  conclusion 
probably  applies  to  any  two-dimensional  transport  problem  requiring 
the  use  of  a  P^  scattering  representation. 
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The  second  basic  problem  associated  with  the  EPFBE  involves  the  odd- 
parity  fluence  transformation.  As  demonstrated  in  Chapter  III  this  trans¬ 
formation  forces  a  certain  anisotropy  on  the  solution  of  a  transport 
problem  and  adversely  affects  the  collocation  evaluation  of  any  term 
in  the  EPFBE  or  boundary  condition  that  contains  the  odd-parity  fluence. 

The  form  of  this  anisotropy  is  fixed  and  at  best  can  only  be  scaled  in 
magnitude  by  the  solution  coefficients  determined  from  the  collocation 
method.  The  only  possible  situation  that  would  allow  an  acceptable, solu¬ 
tion  is  when  the  anisotropy  generated  from  the  odd-parity  fluence  trans¬ 
formation  matches  the  anisotropy  of  the  true  solution.  Such  a  situation 
is  highly  improbable  and  certainly  was  not  the  case  for  the  air-over¬ 
ground  problem. 

Attempting  to  derive  a  synthesis  trial  solution  that  accurately 
represented  the  anisotropy  of  a  particular  problem  would  be  extremely 
difficult.  Such  a  synthesis  function  would  need  to  satisfy  the  constraints 
imposed  by  the  even-parity  fluence  and  be  transformable  using  derivative 
and  integral  operations  into  a  function  that  accurately  approximates  the 
odd-parity  fluence  or  anisotropy.  The  central  difficulty  of  this  deriva¬ 
tion  would  be  understanding  what  constitutes  a  good  representation  of  the 
anisotropy  for  a  certain  problem.  Most  transport  problems  are  so  complex 
that  accurately  approximating  the  anisotropy  would  be  virtually  impossible. 
Additionally,  the  type  of  transport  problem  that  the  above  derivation 
would  be  applicable  must  have  a  fixed  anisotropy  across  the  spatial  domain. 
Few  realistic  problems  meet  this  constraint. 

The  forced  anisotropy  originating  from  the  odd-parity  fluence 
transformation  also  affects  a  pure  spline  trial  solution.  The  effect  of 


this  anisotropy  appears  between  the  nodes  and  modifies  the  shape  of  the 
spline  functions  in  .the  angular  domain.  At  the  nodes  the  collocation 
method  forces  the  correct  solution.  This  modification  from  the  original 
spline  shape  is  the  forced  anisotropy  induced  by  the  odd-parity  fluence 
transformation.  The  only  cure  for  this  problem  is  to  space  the  nodes 
closer  together  in  the  angular  domain.  The  consequence  of  this  action 
is  to  create  a  global  matrix  of  high  order  that  is  costly  to  compute 

and  inefficient  to  solve.  The  optimum  spacing  of  the  nodes  would  have 

*  + 

to  be  determined  through  experimentation  and  would  vary  from  one  type 
of  transport  problem  to  another. 

It  is  concluded  from  this  research  that  the  collocation  method  can 
not  be  effectively  applied  to  the  EPFBE.  This  statement  must  be  qualified. 
The  synthesis  method  is  totally  incompatible  with  the  requirements  of 
forming  the  odd-parity  fluence  through  the  odd-parity  fluence  transformation 
and  thus  an  efficient  solution  method  based  on  synthesis  and  collocation 
is  not  possible.  A  pure  spline  trial  solution  might  generate  an  accept¬ 
able  transport  solution  with  the  collocation  method;  however,  the  amount 
of  discretization  required  in  the  angular  domain  to  offset  the  forced 
anisotropy  would  prohibit  the  classification  of  this  approach  as  efficient. 

The  multigroup  problems  described  in  Chapter  IV  are  common  to  any 
solution  technique  proposed  for  use  with  the  EPFBE.  The  structure  this 
method  imposes  compounds  the  complexity  already  inherent  in  the  EPFBE. 

This  complexity  increases  with  the  number  of  energy  groups  and  would 
restrict  the  use  of  the  many  grouped  cross  section  sets  commonly  used 
in  an  air-over-ground  problem.  Accurate  transport  solutions  are  not 
possible  without  using  these  cross  sections.  A  conclusion  of  this 
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research  is  that  the  standard  multigroup  method  can  not  be  used  with  the 
EPFBE  formulated  for. an  air-over-ground  problem. 

Three  significant  points  concerning  the  EPFBE  have  resulted  from 
this  research.  Miller,  Lewis,  and  Rossow  (Refs.  25  and  26)  presented 
a  solution  technique  that  successfully  eliminated  the  ray  effect.  The 
finite  element  method  applied  to  the  even-parity  functional  appeared  to 
have  potential  for  other  types  of  transport  problems  that  are  affected 
by  this  numerical  problem.  This  research  has  definitively  demonstrated 
that  the  work  of  the  above  authors  can  not  be  extended  to  an  air-over¬ 
ground  problem  and  is  probably  limited  to  the  simple  geometries  and  scat¬ 
tering  models  used  in  their  published  work. 

This  research  has  also  provided  the  basic  criteria  for  attempting 
future  work  with  the  EPFBE.  The  solution  technique  used  to  solve  this 
equation  must  be  integrally  defined  and  contain  the  vacuum  boundary 
condition  as  a  natural  condition.  The  trial  solution  selected  for  use 
with  the  EPFBE  must  be  globally  defined,  thus  eliminating  the  inefficient 
numerical  integration  problems  encountered  in  Chapter  II.  If  the  trial 
solution  contains  an  angular  synthesis  function  it  must  be  simple  in  form 
and  allow  those  parameters  related  to  the  pitch  of  the  streaming  direc¬ 
tion  to  be  separated  from  other  angular  variables. 

The  third  point  addresses  the  problem  of  reconstructing  the  Boltzmann 
fluence  from  a  numerical  solution  of  the  EPFBE.  As  demonstrated  in  Chapter 
III,  the  odd-parity  fluence  transformation  does  not  accurately  represent 
the  anisotropic  portion  of  the  Boltzmann  fluence.  The  central  problem 
with  this  transformation  is  its  dependence  on  the  spatial  and  angular 
derivatives  of  the  selected  trial  solution.  Typical  finite  element  trial 
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solutions  are  selected  primarily  on  the  requirement  of  mathematical  con¬ 
tinuity.  Synthesis  trial  solutions  are  selected  because  they  accurately 
approximate  some  functional  variation  of  interest.  For  the  even-parity 
equation  neither  of  these  types  of  trial  solutions  are  satisfactory 
because  their  angular  and  spatial  derivatives  do  not  accurately  approxi¬ 
mate  the  anisotropic  portion  of  the  Boltzmann  fluence.  Additionally, 
the  specification  of  the  EPFBE  does  nothing  to  constrain  these  deriva¬ 
tives  for  accurately  approximating  the  odd-parity  fluence.  It  is  a  con- 

*  + 

elusion  of  this  research  that  the  anisotropic  component  of  the  Boltzmann 
fluence  can  not  be  accurately  calculated  using  standard  finite  element 
and  synthesis  techniques.  The  impact  of  this  conclusion  is  that  the 
Boltzmann  fluence  can  not  be  accurately  derived  from  a  solution  of  the 
EPFBE  resulting  from  the  use  of  these  numerical  methods.  Successful 
work  using  this  equation  was  always  reported  in  terms  of  the  scalar 
fluence.  The  angular  fluence  results  were  never  presented. 

Recommendations 

Three  facets  of  future  work  are  proposed  for  calculating  scalar 
fluences  with  the  monoenergetic  form  of  the  EPFBE.  First,  a  parametric 
study  should  be  performed  to  determine  the  effect  of  reducing  the  order 
of  the  Legendre  expansion  representing  the  differential  scattering  cross 
section.  If  a  P-j  representation  could  produce  acceptable  results,  the 
decrease  in  complexity  would  allow  the  analytical  formulation  of  the 
weak  form  of  the  EPFBE  using  a  pure  linear  Lagrange  trial  solution. 

Second,  the  angular  synthesis  function  presented  in  Chapter  III  for 
use  with  the  collocation  method  should  be  attempted  with  the  weak  form 
of  the  EPFBE.  The  successful  analytical  formulation  of  the  EPFBE  using 
the  collocation  method  grants  much  credence  to  this  proposal.  Additionally 
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this  angular  synthesis  function  allows  the  pitch  parameter  to  be  separated 
from  the  other  angular  variables.  This  property  would  permit  a  total 
separation  of  the  spatially  and  angularly  dependent  terms  and  thus  the 
individual  integral  evaluation  of  each.  Since  this  angular  synthesis 
function  is  globally  defined  all  analytical  evaluations  could  be  accom¬ 
plished  using  definite  integrals.  These  definite  integrals  would  not 
create  the  large  number  of  terms  in  MACSYMA  experienced  in  the  indefinite 
integration  process  used  in  Chapter  II.  This  proposal  is  not  without 
risk.  The  trigonometric  functions  in  the  angular  synthesis  function 
must  be  multiplied  by  similar  functions  as  specified  by  the  Galerkin 
method.  These  products  of  trigonometric  functions  are  analytically 
evaluated  using  reduction  formulas  that  can  generate  several  terms. 

This  increased  number  of  terms  coupled  with  the  extra  integration  over 
the  angular  domain  required  by  the  Galerkin  method  could  cause  problems 
with  MACSYMA’s  limited  memory.  If  memory  problems  evolve,  a  combination 
of  analytical  and  numerical  integration  could  prove  feasible. 

The  third  proposal  involves  a  modified  form  of  the  synthesis  trial 
solution  presented  in  Chapter  III.  This  trial  solution  can  be  modified 
for  use  with  the  first-order  Boltzmann  equation.  Adding  an  odd  function 
to  this  trial  solution  would  result  in 
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and  all  other  terms  have  been  previously  defined.  This  odd  term  allows 
the  anisotropy  of  the,  first-order  Boltzmann  fluence  to  be  defined  and 
specifically  controlled  by  the  coefficient.  This  control  was  not 

available  in  the  EPFBE. 

The  trial  solution  in  Eq.  (5.1)  could  be  used  with  several  formula¬ 
tions  of  the  first-order  Boltzmann  equation.  The  collocation  method 
could  be  directly  applied  to  this  equation,  since  no  odd-parity  fluence 
transformation  exists.  The  one  possible  drawback  to  applying  this  method 
involves  the  hyperbolic  nature  of  the  first-order  Boltzmann  equation. 
Articles  do  not  appear  in  the  literature  that  successfully  demonstrate  the 
application  of  the  collocation  method  to  a  hyperbolic  equation.  The 
work  by  Houstis,  Lynch,  Rice,  and  Papatheodorou  (Ref.  35)  considered 
only  elliptic  equations.  A  research  project  attempting  this  application 
should  initially  study  the  effects  of  collocating  hyperbolic  equations. 

Another  numerical  method  that  would  be  compatible  with  the  trial 
solution  in  Eq.  (5.1)  and  first-order  Boltzmann  equation  is  the  weighted 
residual  technique  used  by  Roberds  and  Bridgman.  Their  trial  solution 
had  the  same  number  of  unknowns  per  spatial  point  as  the  one  in  Eq.  (5.1), 
but  only  provided  a  P-j  type  of  angular  approximation.  The  ///  term  in 
Eq.  (5.1)  offers  a  P^  type  of  approximation  for  the  same  number  of 
unknowns.  Roberds  and  Bridgman  used  the  method  of  finite  differences 
to  evaluate  the  spatial  varition  of  the  fluence.  The  splines  in  Lq. 

(5.1)  would  serve  the  same  purpose  and  could  be  analytically  evaluated 
using  MACSYMA. 


82 


BIBLIOGRAPHY 


1.  Bell,  G.I.  and  Glasstone,  S.,  Nuclear  Reactor  Theory,  New  York:  Van 
Nostrand  Reinhold  Co.,  1970. 

2.  Straker,  E.A.,  The  Effect  of  the  Ground  on  the  Steady-State  and  Time- 
Dependent  Transport  of  Neutrons  and  Secondary  Gamma  Rays  in  the 
Atmosphere,  Nucl ,  Sci.  Eng.,  46:344-355  (1971). 

3.  Shulstad,  R.A.,  An  Evaluation  of  Mass  Integral  Scaling  as  Applied 
to  the  Atmospheric  Radiation  Transport  Problem,  DS/PH/76-3,  Wright- 
Patterson  AFB,  OH:  AFIT,  November  1976. 

4.  Mynatt,  F.R.,  Muckenthal er ,  F.J.,  and  Stevens,  P.N.,  "Development 
of  Two-Dimensional  Discrete  Ordinates  Transport  Theory  for  Radiation 
Shielding,"  Union  Carbide  Corporation,  Nuclear  Division  (1969)*. 

5.  Carlson,  B.G.  and  Lathrop,  K.D.,  "Transport  Theory  -  The  Method  of 
Discrete  Ordinates."  Computing  Methods  in  Reactor  Physics,  Chapter 
III,  Greenspan,  Kelber,  and  Okrent,  Eds.,  Gordon  and  Breach,  New 
York  (1968). 

6.  Lathrop,  K.D.,  "Ray  Effects  in  Discrete  Ordinates  Equations."  Nucl . 

Sci,  Eng.,  32:357-369  (1968). 

7.  Lathrop,  K.D.,  "Remedies  for  Ray  Effects,"  Nucl .  Sci  ,  Eng . ,  45:255- 
268  (1971). 

8.  Jung,  Jungchung,  "Discrete  Ordinate  Neutron  Transport  Equation  Equi¬ 
valent  to  PL  Approximation,"  Nucl .  Sci .  Eng. ,  49:1-9  (1  972). 

9.  Reed,  W.H.,  "Spherical  Harmonic  Solutions  of  the  Neutron  Transport 
Equation  from  Discrete  Ordinate  Codes."  Nucl.  Sci.  Eng.,  49:10-19 
(1972). 

10.  Jung,  Jungchung,  "Second  Order  Discrete  Ordinate  P.  Equations  in 
Multi-Dimensional  Geometry."  Journal  of  Nuclear  Energy,  27:577-590 
(1973). 

11.  Jung,  Jungchung,  "Solution  of  Standard  Diamond  Difference  Equations 
for  Discrete-Ordinate  Neutron  Transport  Equations  Equivalent  to  the 
PL  Approxmation  in  x-y  Geometry."  Nucl .  Sci .  Eng . ,  53:355-364  (1974). 

12.  Lathrop,  K.D.,  Private  Communication.  Energy  Division,  Los  Alamos 
Scientific  Laboratory,  January  17,  1978. 

13.  French,  R.L.,  "A  First-Last  Collision  Model  of  the  Air/Ground  Inter¬ 
face  Effects  on  Fast  Neutron  Distributions,"  Nucl  ,  Sci .  Eng. ,  19:151- 
157  (1964). 

14.  French,  R.L.,  Mooney,  L.G.,  "Last-Collision  Model  of  the  Air-Ground 
Interface  Effect  on  Fast-Neutron  Angle  Distributions,"  Nucl .  Sci .  Eng. , 
19:1  51-157  (1  964). 


83 


15.  Guy,  R.W.,  Albert,  T.E.,  and  Straker,  E.A.,  "A  Comprehensive  Model 
for  Initial  Weapons  Radiation  Environments,"  Tran.  Am.  Nucl .  Soc . , 

1  9:452  (1  974). 

16.  Vladimirov,  V.S.,  "Mathematical  Problems  in  the  One-Velocity  Theory 
of  Particle  Transport,"  Atomic  Energy  of  Canada,  Ltd.,  Ontario  (1963), 
translated  from  Trans.  V.A.  Steklov  Mathematical  Institute,  61  (1961). 

17.  Kaplan,  S.  and  Davis,  J.A.,  "Canonical  and  Involuntary  Transformation 
of  the  Variational  Problems  of  Transport  Theory."  Nucl.  Sc i .  Enq.,  28: 
166-176  (1967). 

18.  Kaplan,  S.,  Davis,  J.A.,  and  Natelson,  M.,  "Angle-Space  Synthesis--An 
Approach  to  Transport  Approximations."  Nucl.  Sci .  Enq.,  28:364-375 
(1967). 

19.  Natelson,  M.,  "A  Strategy  for  the  Application  of  Space-Angle  Synthesis- 
to  Practical  Problems  in  Neutron  Transport."  Nucl.  Sci.  Enq.,  31:325- 
336  (1968). 

20.  Kaplan,  S.,  "Synthesis  Methods  in  Reactor  Analysis."  Advances  in 
Nuclear  Science  and  Technology,  3:233-266  (1966). 

21.  Stacey,  W.M.,  Jr.,  "Flux  Synthesis  Methods  in  Reactor  Physics." 

Reactor  Technology,  15:210-238  (Fall  1972). 

22.  Kaplan,  S.,  "Some  New  Methods  of  Flux  Synthesis."  Nucl ,  Sci .  Eng . , 
13:22-31  (1962). 

23.  Roberds,  R.  and  Bridgman,  C.,  "Application  of  Space-Angle  Synthesis 

to  Two-Dimensional  Neutral-Particle  Transport  Problems  of  Weapon 
Physics."  Nucl  .  Sci.  Enq.,  64:332-343  (1977). 

24.  Miller,  W.F.,  Lewis,  E.E.,  and  Rossow,  E.C.,  "The  Application  of 
Phase-Space  Finite  Elements  to  the  One-Dimensional  Neutron  Transport 
Equation."  Nucl .  Sci .  Eng. ,  51:148-1  56  (1973). 

25.  Miller,  W.F.,  Lewis,  E.E.,  and  Rossow,  E.C.,  "The  Application  of 
Phase-Space  Finite  Elements  to  the  Two-Dimensional  Neutron  Transport 
Equation  in  x-y  Geometry."  Nucl  ,  Sci  .  Enq . ,  52:12-22  (1973). 

26.  Reed,  W.H.,  Hill,  R.R.,  Brinkley,  F.W.,  and  Lathrop,  K.D.,  "TRIPLET: 

A  Two-Dimensional,  Multigroup  Triangular  Mesh,  Planar  Geometry, 

Explicit  Transport  Code."  Los  Alamos  Scientific  Laboratory,  LA- 5428- 
MS,  October  1973. 

27.  Kaper,  H.G.,  Leaf,  G.K.,  and  Lindeman,  A.J.,  Applications  of  Finite 
Element  Methods  in  Reactor  Mathematics.  Numerical  Solution  of  the 
Neutron  Transport  Equation.  ANL-8126,  Argonne  National  Laboratory 

(1974T. 

28.  Briggs,  L.L.,  Miller,  W.F.,  and  Lewis,  E.E.,  "Ray-Effect  Mitigation 
in  Discrete  Ordinate-Like  Angular  Finite  Element  Approximations  in 
Neutron  Transport,"  Nucl .  Sci .  Eng. ,  57:205-217  (1975). 


84 


29.  Zachmanoglou,  E.C.  and  Thoe,  D.W.,  Introduction  to  Partial  Differential 
Equations  with  Applications,  Baltimore:  The  Williams  and  Wilkens 
Company,  1976. 

30.  Luenberqer,  D.G.,  Optimization  by  Vector  Space  Methods,  John  Wiley 
Co.,  1969. 

31.  Huebner,  Kenneth  J.,  The  Finite  Element  for  Engineers,  John  Wiley 

Co.,  1  975.  "  '  '  '  “ 

32.  Wheaton,  Ronald  C.,  AFIT/GNE/PH/78D-1 5,  Wright-Patterson  AFB ,  OH: 

AFIT,  December  1978. 

33.  Strang,  Gilbert,  and  Fix,  George  J.,  An  Analysis  of  the  Finite  Element 
Method,  Prentice-Hall  Inc.,  1973. 

34.  MACSYMA  Reference  Manual,  Massachusetts  Institute  of  Technology,  1977. 

35.  Houstis,  E.N.,  Lynch,  R.E.,  Rice,  J.R.,  and  Papatheodorou ,  T.S., 
"Evaluation  of  Numerical  Methods  for  Elliptic  Partial  Differential 
Equations,"  Journal  of  Computational  Physics,  27:323-350  (1978). 

36.  Davis,  Harry  F.,  Introduction  to  Vector  Analysis,  Allyn  and  Bacon, 

Inc.,  1967. 


* 


APPENDIX  A 


DERIVATION  OF  EVEN-PARITY  FORM  OF  THE 
BOLTZMANN  NEUTRON  TRANSPORT  EQUATION 
The  monoenergetic  steady-state  Bo-ltzmann  neutron  transport  equation 
can  be  written  as 

A.  c^a)  i-c7r(ft  J -  J  j  ctt  - Sct.Aj  -  O 

Jt*n 

Since  this  equation  holds  for  A.  it  must  also  hold  for -A 

-A-vlUc?-*) -fazC*,*-*)  Jc^Vo^A-  Sc?,.*)  =  O'*  (A. 2) 


If  Eqs.  (A.l)  and  (A. 2)  are  added  together  the  following  results: 

A.  v(FiZA)  -  J(v aW  +  oji*)  (£(?,».)*■  jEi*,-**)  ~ 

j  <*)*'••  Aj ) SlWdK-  J  =  O 

jhh 

Substracting  Eq.  (A. 2)  from  Eq.  (A.l)  results  in 

A>  V  (  £ir,/L)  +1^-*.))  +oy{r/  (1V,aJ- J i?r\})  ~ 

J  (  C£  a.'\J  -  -A  I)3lC',  A.'Jcl%-  tScPA.1  ~  )  =  O 

'L.  w 


(A. 3) 


(A. 4) 


Now  define  the  following  quanties: 

H^aj  =  2  LJcf;a/ «■  Jcr-4/] 


(1.10) 


3Ccf;A)  3  2  ChwIcAA'J 


(A. 5) 


ov^aa;*  2  Ccj^a  am  c^crA-'xJ  1 


(1.16) 
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r 

Cs^cp'Si'A.)  =  2  cr(rA'-sijJ 


S^l^A.)  -  2  [  Scf^/iJ-  SjlP,-A)~] 

53ik*j  =  i  [  Sws<K-*,J 

Substituting  these  definitions  into  Eqs.  (A. 3)  and  (A. 4)  respectively 
results  in 

A-  V  XCr.Al  +  0>(r;  Vi?,*)-  j  O^C^a'-T;  Xo~A'j  c/a' 

-/V/t  J 

-  )  =  O 

A  v  ♦  aTct)  X( K*J  -  I <%u(?, Z'ZjTcr, X)  d«.‘ 

An 

~  Cf^A)  =  O 


(A. 6) 


(1.11) 


(1.12) 


(A. 7) 


(A. 8) 


Now  define  the  following  operators 


-fc?A)  -  <7r<r;  fc^/U  -  (  O^W-Ai  'fcfi,*)  dX 

An 

Gm  -fcr  A)  -  (rr(r)  fee,*.)- f  -f  <.?,*)  U*‘ 


(1.14) 


(A. 9) 


To  use  the  above  two  equations  the  following  must  be  proved 

{  Cu,  £ir,ZjdA.'  =  { ozur,  *'<J  Vl?  x  )cU'  (A. 10) 

Hn  J  An 

This  can  be  done  by  using  the  previously  presented  definitions 

j  'i  ) )  (iiPA  i j )  d*. '  (A. 11) 

JHn  J 
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The  term  on  the  right  hand  side  of  Eq.  (A. 11)  can  be  expanded  to  give 


*1  /  (  o-5(r*  A'AI  It w  *  q  c r,A'-A/  !(?,*)  *  Itf-XJ  -* 

'sn 

+ozC?,A‘-/ilIlK-Z))  d*  (A. 12) 


however,  the  following  relationships  are  true 
T  j  c$c?A-;uIci-*)cti'  =  ^  f  GCKf~*J  IcK-X)  clA‘ 

An  An- 


‘if  Hir.Xtdf  *  H  j  c/A‘ 

JHn  An 


(A. 13) 


(A. 14) 


Substituting  the  results  of  Eqs.  (A;13)  and  (A. 14)  into  Eq.  (A. 12)  gives 

7/  (  qfr'  a'A;  +q^A‘-/l» )  Hc?,A‘l  UA-  (A. 15) 

An 


Using  Eqs.  (1.16)  and  (A. 15)  gives 

/  OIoC^  a'  a;  dA' 

/47T 

thus  proving  the  equality  in  Eq.  (A. 10).  A  similar  proof  can  be  con¬ 
structed  to  show 

(  01u,l  ?,£•*)  Ut'  =  /  Oj^CF'A'AJ  XirA'iUA'  (A. 16) 

An  An 


Using  the  results  of  Eqs.  (A. 10)  and  (A. 16)  with  the  definitions  of 
Eqs.  (1.14)  and  (A. 9)  in  Eqs.  (A. 7)  and  (A. 8)  gives 

Gcj  =  >$3  (■?.*)  -  a- V  Xcr,At 


(A. 17) 


GmX^aI*  rt v  fcKAJ 


(A. 18) 


Eqs.  (A. 17)  and  (A. 18)  represent  a  first  order  set  of  coupled  integro- 
differential  equations.  Now  define  two  new  operators  such  that 
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(A. 19) 


K<j  Gj  f(^A)  =  -ft*,*) 

K^GU  -  ft?. a>- 

Kg  and  Ku  are  the  inverse  operators  of  Gg  and  Gu  respectively.  The 
derivation  of  these  two  operators  is  presented  in  Appendix  B.  Applying 
these  inverse  operators  to  Eqs.  (A. 17)  and  (A. 18)  results  in 


Vc^a)  = 

Kj  (  SjCf.A  l  -  A-  V  XCr,Al  ) 

(A. 21) 

yi*.  a)  = 

KA  ( -  A  v  Vi?,*) ) 

(A. 22) 

Substituting  Eqs.  (A. 21)  and  (A. 22)  into  Eqs.  (A. 17)  and  (A. 18)  gives 
Av  KPLSiL(ftA/-A  vkU'(A'vfttA>)  +  G3i'<':*>-S9(tA>:  O  (a.23) 

A-  v  -  A-v  Kj  (Av  X(?,a>)  «■  S^v-  O  (A. 24) 

Eq.  (A.23)  is  called  the  second  order  even-parity  form  of  the  Boltzmann 
neutron  transport  equation.  Eq.  (A. 24)  is  the  odd-parity  form  of  this 
same  equation.  Rearranging  Eq.  (A.23)  gives 

-/V  V  ku  (A- V  Vet, A) )  *  Gj  Ycfi.i./  =  A  v  k^S^Cr.A.)  (1.9) 
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APPENDIX  B 


-DERIVATION  OF  THE  Ku  OPERATOR 


The  Ku  operator  as  seen  from  Appendix  A  takes  the  inverse  of  the 
Gu  operator.  The  Gu  operator  is  defined  as 

Qu.  ~PcAi  -  crT-f(A.)  -  I  A’K )  -f(A')  clA' 


(B.l ) 


(The  t  variable  has  been  omitted,  since  it  does  not  effect  this  deriva¬ 
tion.)  Let  be  defined  in  terms  of  a  Legendre  polynomial 


expansion 


jL*  o 


(B.2) 


From  the  addition  theorem  of  Legendre  polynomials  the  following  can  be 


stated 


*  2£+  ,  y/n  (A'J 


(B.3) 


where  Yjtf.)  is  a  normalized  spherical  harmonic  and  Y^xV)  is  the  complex 

conjugate  of  \lA)  .  Substituting  Eq.  (B.3)  into  Eq.  (B.2)  yields 

oo  X  t 

Ca'aI  =  £Z  Z£  OjM  (B 


Now  assume  that  ~F(A)  can  be  represented  in  a  spherical  harmoni 


OO  JL 


-Fu)  =  Z  Z  £  Yt  u  i 


(B.4) 


c  expansion 
(B.5) 


Placing  Eqs.  (B.4)  and  (B.5)  in  Eq.  (B.l)  gives 

r  •  4  /■„  i  oo  J 

G  TCAI*  OlZZ"  U)  -  ZI  C^Y.^YaKiTT  -f;  YLkIA‘J  do.*  (B.6) 


Rearranging  Eq.  (B.6)  leads  to 


.  i.  i"  mmm 


Gu-f(  A  ):o-TIZ  fx  )-ZZq£  2  f  Y  CaV  Xcxi  cl*.'  (B .  7 ) 

Tj ■*'  .<*»  «-.-4  ■*  yjrOit..j  J  >A*  ■*«. 


The  orthogonality  properity  of  spherical  harmonics  states 

f"  )iU  )  \,k  ( A  ^  ~  « 

S7J 


(B.8) 


where 


4, 


■  1  'f 


’J  {  --  O  ;f  it. 


Using  the  orthogonality  property  in  Eq.  (B.7)  results  in 

Gvfc,= 


(B.9) 


Combining  terms  in  Eq.  (B.9)  gives 

G-fw,=  tl-L 


(B. 10) 


Eq.  (B.10)  illustrates  that  the  Gu  operator  alters  the  spherical  harmonic 
coefficient  ^  by  a  factor  of  «rT- <z«J  .  By  definition  an  inverse 


operator  to  must  do  the  following 


Gu  -f(A )  T  f(Aj 


(B.ll) 


Therefore  if  Ku  is  defined  as 


(o-r  -  '<Z«J 


(B.12) 


then 


'*■  ^  A* O  *■ 


thus  satisfying  Eq„  (B.ll) 


91 


Eq.  (B.12)  can  be  rearranged  into  the  following  form: 

^  mlL (A  *  A  ( -rfr-)  *  W  •  lB-’4) 

Rearranging  Eq.  (B. 14)  and  again  using  the  orthogonality  property  of 
spherical  harmonics,  the  following  results: 

Kjc*)  =  oy  a  I  rjJ.  £.j  fj  (B'15) 


Now  make  the  foil 


Oiu  M 


■■V  II 


owing  definition 
op  _±  crf- 


Jt=0  **  =  •£. 


(B.16) 


Using  the  addition  theorem  on  Eq.  (B.16)  gives 

oo  /  2JL+  i  )  (B.17) 

<WavC)  =  2.  oy-$~  t  /  If  (**/ 

1=0 

Using  the  results  of  Eq.  (B.16)  and  previous  definitions  the  following 
can  be  stated: 

Ku  -fui  -  oy  (  r 

Using  a  similar  derivation  the  following  form  of  Kg  can  be  determined: 

k3  -fo)  =  cy  (  *  jT  ^3  -fw  cU’  )  (B.  18) 


o-Ku(Wf(, 


'iclK') 


(1.13) 
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APPENDIX  C 


DERIVATION  OF  THE  WEAK  FORM  OF 
SECOND  ORDER  EVEN-PARITY  FLUENCE  EQUATION 
The  classical  technique  for  deriving  the  weak  form  of  an  equation 
can  be  illustrated  in  the  following  general  example.  Let  the  equation 
be  represented  by 


L  Ycx)  =  -fc?; 


(C.l) 


where  L  is  an  operator,  Yet)  is  the  dependent  variable,  -fctj  is  the 
source  term,  and  "x  is  a  vector  representing  the  independent  variables. 

If  Eq.  (C.l)  is  rearranged  and  multiplied  through  by  a  weight  function, 

CCij  ,  the  following  results: 

(  L  Yet)  ~  fct)  )  ZCx)  *  O  ( c *2) 


Integrating  Eq.  (C.2)  over  its  domain  of  definition  gives 

f  (LYcxJ-  ~f(xJ  )  Ux)  dx  =  O  (C-3) 

‘'o? 

Eq.  (C.3)  is  the  weak  form  of  Eq.  (C.l).  This  equation  is  "weak"  in 
that  it  holds  in  an  integral  sense  rather  than  a  point-wise  sense  as  in 
Eq.  (C.l). 

The  derivation  of  the  weak  form  of  the  even-parity  equation  follows 
the  above  illustration.  The  even-parity  equation  is  given  by 
-a.-v  A  v  *  Gj  A  (c.4) 


where  the  operator  for  this  equation  is 

L  =  -A.v  Ku  A-v  *  G3 


(C.5) 
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and  the  source  term  is 


■f CP, A)  =  5*X T,ai  -  A- o  Ku  (^4  I 


(C.6) 


Rearranging  Eq.  (C.4)  and  multiplying  through  by  the  weight  function 
results  in 

(-  A- v  Ku  /i-v  *  Gj  -  SjcK'Cj  *A v  )  =  O  (C.7) 


Now  integrate  over  the  phase  space  defined  by  "r  and  .A 
[(-*.  v  Ku  a.  v  fcf'A  j)  £(?,*.  I  *  CGy  )  ec?,*.!  - 

S^a.)  Ci?,A.)  ♦  (*:y  Ku5^t?,Zj)  ciZd?  -  O 


(C  .8) 


For  simplification  the  following  definition  is  presented 
<  Juft),  *  j  A ctj.)  ct* 

''L 

Substituting  Eq.  (C.9)  into  Eq.  (C.8)  results  in 

C  <-/vv  a.  v  £C?,AI?  t  <G cjVcr.A.l.CCt'A.i)  ~ 

<  ec?,Z)>  +  <A^Ku5mc^i,  £ cfi,A.,>]dr  -  O 

Again  for  simplification  let 

C^C  ^  A  )  =  KH  A-  V  ft?, A  J 


(C.9) 


(C.10) 


(C.ll) 


and  operate  on  each  term  in  Eq.  (C.10)  individually.  Using  Eq.  (C.ll) 
in  the  first  term  of  Eq.  (C.10)  results  in 

~  <  A- V  <£(?,*),  S <.?,*.)>cL?  (C.  12) 

Since  A  and  v  have  independent  coordinate  systems,  Eq.  (C.12)  can  be 
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modified  to 


-j~  <  V.-A^fF.Ai.CCI’Aj  >cLt 


(C. 13) 


Using  the  following  vector  identity  (Ref.  36) 

V  CCj“Aj  (A^C^a.i)  -  £ci?A.I  V ■  (AcfCZ*.))  -  (a vccf.aj 


£Cff/tJ  V  (/Left?,  *>) s  V-  £lP,Al(A<llZto  b  (A  tpP.AI)  •  VU?.*) 


(C. 14) 


(C. 15) 


in  Eq.  (C.13)  results  in 

-  (v- A^cf'A)  simdtcU  =  £fjv- * 

A  ^  \  t.  ( C . 

j'  j' (ac^C?a)'  Vf cL ’?  cJ-A 

A  £ 

Applying  the  divergence  theorem  to  the  first  term  on  the  right  hand  side 


(C.16) 


of  Eq.  (C.16) 


~  j'  £C?aI  <£(5^  a  I  els  cLa. 


(C.17) 


-where  S  represents  the  boundary  of  the  spatial  domain.  Substituting 
the  definition  of  c^aJ  into  Eq.  (C.16)  and  using  the  results  of  the 

divergence  theorem  modifies  the  first  term  of  Eq.  (C.10)  to 

j  K.  A  V  A'Vf^,Al)c(r  —  |  (.A-  n. )  (  Ku  V  'fc.r,  A/  ^  o/-S  (C  .  18  ) 

Operating  on  the  fourth  term  in  Eq.  (C.10)  with  both  the  vector  identity 

in  Eq.  (C.15)  and  the  divergence  theorem  results  in  the  following  new 
form 

-j  <  a/(  /v  oU?A.)>d?  +•  fU,  K^SMi^Ai>df  (C.19) 


The  weak  form  of  Eq.  (C.4)  is  now  given  by 
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(C.20) 


< KuSh(?,x. i,Zva?,A) >Jd? *■ if*'* j, KM i- a* v H{^A.))>ds  -  O 

By  the  definition  of  the  odd-parity  fluence  presented  in  Eq.  (A. 22), 
the  term  in  Eq.  (C.20)  representing  and  integration  around  the  spatial 
domain  can  be  modified  to 

[  <  £c?,  a)  (A-n.),  Xc?,A)  ?  c/s  (C.21) 

J5v 


Assuming  a  vacuum  boundary  condition  the  following  is  true 


W,AJ  + 

XCrtA.  )  '  - 

O  ,  ^  e-S^, 

,  A- A  <  O 

(C.22) 

- 

'XCrlA  1  = 

0  ,  ^eS'v  , 

,  A’  A  >  0 

(C.23) 

Eqs.  (C.22)  and  (C.23)  can  be  rewritten 
Wcr.AI  =  -  XCP,A  I  ,  reSv, 

Vc?z.)  =  X  cr, a.)  ,  r  e  S„  , 


as 

a- a  <  O 

(C.24) 

A- a  >0 

(C.25) 

Using  the  results  of  Eqs.  (C.24)  and  (C.25)  the  following  modification 
to  Eq.  (C.21)  can  be  made: 


j  <  CC?A.)  CA-A),  A,  A )  >  ds  -  j  j-  CCfA.)  %?*.)  (AA)dAcli 

-WV  -T  4  _ 


o 

Eq.  (C.26)  is  equivalent  to  the  following 


f  [  £c?,Al  (A- Velvets' 


(C.26) 


f  <  UKA)  lA-Kl,  V<.?,Ai>ds 
% 


(C.27) 
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Thus  the  weak  form  of  the  even-parity  equation  that  naturally  satisfies 
the  vacuum  boundary  condition  is 


a-^  /vv  CCf.A.)?  +  <  VcZAj,  L:rtA.i  >  ~ 

^  £c?,  a  )  I  A.' n  I J  ^tr,^  J  y  els  —  O 


(2.1) 
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APPENDIX  D 


•BILINEAR  LAGRANGE  POLYNOMIALS 
A  bilinear  Lagrange  polynomial  can  be  represented  as 

Lcx.y)  =  a,xy  +-a;x  +  Q^y 

If  a  rectangular  element  is  defined  as  shown  in  Fig.  16,  Eq.  (D.l)  can 
be  properly  constrained  to  give  the  tent  function  used  as  a  trial  solu¬ 
tion  in  Chapter  II.  Let  Lp  L^,  L3,  L4  be  the  value  of  Eq.  (D.l)  at  the 
nodes  defined  in  Fig.  .  Then  the  following  system  of  equation  can  be 
written 

(D.2) 

L.  ( x.,y.)  *  a,x,y,  tq^x,  t  Qjy,  t  qh 

Lj  =  a.Xjy,  j  t  a., 

Lmcx„yj=  a.vys,  *a,x,  +  +  aw 


Now  define  the  matrix  H  such  the 


x.Y. 

X, 

y* 

1 

x*y  i 

V* 

1 

y* 

Xj 

Y? 

1 

*1 

/ 

and  the  vector  A  and  L  where 


(D.3) 


Fig  16.  Rectangular  Element 


A  ~  C.  & ,t  d i  i  Ctj  («J 


(D.4) 


V  -  LL„l*.4.U] 


The  system  of  equations  represented  in  Eq.  (D.l)  can  now  be  represented 


[h1  a--  l 


(D.5) 


Multiplying  the  inverse  of  the  H  matrix  on  both  sides  of  Eq.  (D.5)  yields 

A  ~  [h!  L  (D.6) 


Eq.  (D.l)  can  be  expressed  as  a  product  of  a  row  vector  P  and  the  column 


vector  a,  if 


PT  =  C  xy,  X  f  y.  J  U 


(D.7) 


Lcx.vl- 


(D.8) 


Substituting  the  values  of  A  from  Eq.  (D.6)  into  Eq.  (D.8)  gives 
LtK.yJ  =  PT[H]T 


(D.9) 


By  letting  the  L  assume  the  following  values  the  tent  function  of  nodes 
1,  2,  3,  4  can  be  found  respectively: 


C  I.O.  0.0  I 
(  O,  I  ,  o  ,  o  ) 
(0,0,  I  .0  ) 

(  O  ,0,  o,  I  J 
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APPENDIX  E 


...  CANONICAL  TECHNIQUES 
CANONICAL  ELEMENTS/INTEGRATION 

The  concept  of  a  canonical  element  is  used  in  the  finite  element 
method  to  simplify  the  integral  evaluations  associated  with  the  chosen 
numerical  technique.  In  Fig.  17  a  standard  rectangular  element  with 
sides  parallel  to  the  coordinate  axis  (x,y)  is  presented. 


Fig.  17.  Standard  Rectangular  Element 
Inscribed  in  this  element  is  a  natural  coordinate  system  defined  by  (n,  ). 
The  relationship  between  the  natural  coordinates  and  axis  coordinates 
is  given  by 


/ 


(E.l) 


0. 


y-  Vc 
b 


(E.2) 


This  coordinate  transformation  defines  a  canonical  element  that  repre¬ 
sents  the  rectangular  element  in  Fig.  18. 
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A 


Fig.  18.  Canonical  Element 

From  Eqs.  (E.l)  and  (E.2)  the  following  can  be  generated 

dl  _l  _ _ 

dx  ~  a.  ’ 


oi  x  =  excLj 


5*  -  4-  =» 


(E.3) 

(E.4) 


Let  f(x,y)  be  a  function  that  is  identically  defined  over  all  finite 
elements  in  the  (x,y)  coordinate  system  shown  in  Fig.  17.  Using  the 
coordinate  transformation  equations  this  function  can  be  redefined  in 
the  natural  coordinate  system  as  -f(af+  xc,  v«  )  .  Since  this  func¬ 

tion  is  identically  defined  in  all  elements  the  transformation  equations 
will  generate  exactly  the  same  functional  representation  in  the  canonical 
element  for  all  the  elements  in  the  (x,y)  coordinate  system.  The  only 
difference  is  the  scaling  factors  a  and  b  that  are  related  to  the  dimen¬ 
sion  of  the  element.  Thus  the  following  equality  holds 


This  relationship  allows  the  integral  evaluation  of  f(x,y)  over  a  finite 
element  to  be  done  once  over  the  canonical  element.  If  a  finer  discreti¬ 
zation  is  necessary  the  canonical  evaluation  need  be  changed  only  by  the 
new  dimension  parameters  a,  b  caused  by  this  mesh  refinement. 
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GALERKIN  WEAK  FORM  EQUATION/ISOTROPIC  SCATTER 


The  global  integrals  in  the  Galerkin  weak  form  equation  resulting 


from  the  Gg  operator  and  using  isotropic  scatter  is 


^77 


Vtr;A‘jG^' 


(E.6) 


If  is  represented  by  the  trial  solution  of  Eq.  (2.20),  Eq.  (E.6) 

gives 


.  (E.7) 


The  bilinear  Lagrange  polynomials  meet  the  requirement  of  being  identi¬ 
cally  defined  over  all  elements,  thus  the  integral  in  Eq.  (E.7)  can  be 
evaluated  over  a  canonical  element  as 

ab 


)  d.\ci» 


-I  -i 


(E.8) 


The  evaluation  of  Eq.  (E.8)  can  be  done  once  and  the  result  used  for  all 
other  angular  elements. 


GALERKIN  WEAK  FORM  EQUATION/ANISOTROPIC  SCATTERING 
The  global  integral  resulting  from  the  Ku  operator  using  a  P1  aniso¬ 
tropic  scattering  representation  is 


Cliri 


~j.(r)  j  C  /JLU.'+  J  V  1-tA.^CoSCX-X  '))*.  Ct*' 


(E.9) 


If  the  even-parity  fluence  is  again  represented  as  a  tensor  product  of 
bilinear  Lagrange  polynomials,  then  the  derivatives  that  result  from 
the  gradient  operator  are  identical  over  each  angular  element.  The 
integrand  of  Eq.  (E.9)  contains  other  functions  besides  these  derivatives. 
These  functions  vary  differently  over  each  angular  element,  thus  elimin¬ 
ating  a  single  canonical  evaluation  that  would  represent  all  these 
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r  - . 

elements.  Eq.  (E.9)  can  be  evaluated  over  the  canonical  element;  how¬ 
ever,  it  must  be  separately  evaluated  for  each  element  defined  in  the 
angular  domain  since  the  location  of  the  element  now  affects  the  value 


of  this  term. 


APPENDIX  F 


TOTAL  TERMS  IN  A  xjUt.Al 

The  MACSYMA  program  did  not  have  enough  available  core  to  expand 
the  first  term  in  the  Galerkin  weak  form  equation  into  all  its  distinct 
terms.  The  reason  this  core  limit  was  reached  can  be  seen  from  the 
following  evaluation.  The  first  term  of  the  Galerkin  weak  form  equation 
can  be  expressed  as 


j  j  A.  v  I  Kfj.  a  'Pc?,  a .)  c/a 


Using  the  definition  of  the  Ku  operator  expands  the  integrand  in  Eq. 
to  give 


(F.l) 

(F.l) 


r  • 

^  :  i  I 


w  l. 


aTCri  (a- v  *  J  O'^AA)  a  d/i'D  U  (F. 


2) 


For  this  evaluation  only  the  following  is  used 

JL-VWA.)  J  °iu.  CA'  A l  A  -V  WZX I  o(A' 


(F.3) 


Applying  the  definition  of  the  -a  ^  operator  in  cylindrical  geometry 
expands  the  integrand  to 


-A.  VfC 


(  (* 

r.AlJ  CTku  (a  Al\ 

‘in 


•J  >  ~  l  d(pV(r,A) ) 


(F.4) 


f)  J V'  o  i  J  £Xa 

If  f(^Aj  is  represented  by  linear  Lagrange  interpolating  functions  in 

a  four-dimensional  phase  space  (f>  *,*(  * >  the  following  general  expression 

can  be  written 


Aj  s  (a.*  (  e  f  )  (  g  *  A  *  ) 


(F.5) 
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Substituting  Eq.  (E.5)  into  the  first  term  in  brackets  in  Eq.  (E.4)  gives 


'J  i-  W*  ‘cosCx'J 


(tj-dzJ  Ce+f/A.)  (y  +  hx)  Cu+Zbp) 


(F.6) 


Eq.  (F.6)  contains  sixteen  distinct  terms.  Proceeding  in  the  same  man¬ 
ner  with  the  other  two  terms  in  the  brackets  results  in 


C 


a  + 


bp )  (c  +  olt )  ( e  +  -fjx )  ( y  cos  lx  J  +  hsmixJ  +  hys/rilx/) 
JU.  L  ext  bp )  (  e  +  -f ^ )  (.  j* b-x  )  cL 


(F.7) 

(F.8) 


Eq.  (F.7)  contains  twenty-four  distinct  terms  and  Eq.  (F.8)  contains 
eight.  The  total  number  of  distinct  terms  resulting  from  the  bracketed 
terms  in  the  integrand  of  Eq.  (F.4  )  is  forty-eight. 

The  odd-parity  differential  scattering  cross  section  can  be  repre¬ 
sented  for  a  P3  legendre  expansion  as 

crK^tut)  =  01,  “4  *  U?  ( 


where 


M4  =  lxu 


\J t-pL*  J  l-y}  * co  s(*l  Cosin')  *J  1-  pi1'  si'ilxl  Si*(x‘  J  (F.  1 0) 


Eq.  (F.10)  represents  three  distinct  terms  and  cubing  this  equation 
results  in  ten  for  a  total  of  thirteen  distinct  terms.  Therefore,  the 
integrand  in  Eq.  (F.4)  contains  624  distinct  terms. 

For  now  ignore  the  integration  over  A!  and  define  the  weight  function 
as 

CC^AI  =  (aa.  +  bbp)  ice  fdicLi)  (ec  «-  t jy  bkx  )  (F.ll) 
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The  coefficients  are  different  than  those  in  Eq.  (F.5),  because  in  the 
Galerkin  method  different  combinations  of  Lagrange  polynomials  all 
defined  over  the  same  finite  element  are  multiplied  together  (see  Eq. 
(2.18)).  The  A  v c i?,a  )  has  forty-eight  distinct  terms  just  as  a v  / . 
Multiplying  the  forty-eight  distinct  terms  of  Eq.  (F.ll)  by  the  624 
terms  resuHing  from  the  previous  work  generates  a  total  of  29,952 
distinct  terms.  This  number  is  an  absolute  minimum  since  the  integra¬ 
tion  operation  was  not  included.  The  integral  evaluation  over  a  rectan- 

* 

gular  element  of  arbitrary  dimension  would  generate  at  least  four  more 
terms  for  each  term  evaluated.  If  the  integral  was  complex,  then  more 
than  four  terms  could  result. 

Obviously  many  of  the  terms  are  similar  and  could  be  combined  to 
reduce  the  total.  MACSYMA  can  only  begin  this  simplification  after  all 
the  distinct  terms  are  generated.  Several  approaches  were  attempted  to 
alleviate  this  shortcoming  in  MACSYMA,  but  all  proved  unsuccessful.  The 
main  problem  was  the  limited  core  capacity  of  the  DEC  10  on  which  this 
program  resided. 


i 
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APPENDIX  G 


DERIVATION  OF  ANGULAR  SYNTHESIS  FUNCTION 

In  Chapter  III  an  angular  synthesis  function  was  derived  for  use 
with  the  collocation  method.  The  selection  of  this  function  was  guided 
by  two  basic  criteria.  The  first  criterion  dictated  that  the  form  of 
an  angular  synthesis  function  must  be  simple  to  avoid  the  problems 
encountered  in  Chapter  II.  The  second  criterion  is  the  basic  tenet  of 
using  a  synthesis  function  that  accurately  approximates  the  angular 
variation  of  the  neutron  fluence.  As  mentioned  in  Chapter  II,  the 
choice  of  an  appropriate  synthesis  function  is  mai  e  from  experience, 
previous  calculations,  or  intuition. 

The  collocation  angular  synthesis  function  was  selected  from  the 
definition  of  the  even-parity  fluence.  This  fluence  can  be  represented 
as 

i'CrA.I-  2  [  let*/ *  Ict-Ail  (G.l) 

The  work  accomplished  by  Roberds  and  Bridgman  (Ref.  23)  demonstrated 
that  a  spheroid  can  be  used  as  an  angular  synthesis  function  for  the 
regular  Boltzmann  fluence.  Fig.  19a  illustrates  a  spheroid  that  could 
be  used  to  represent  JEcr.A.j  for  some  anisotropic  angular  distribution. 

Fig.  19b  is  a  polar  plot  of  Tc~-*l  using  the  same  spheroid  synthesis 
function.  If  the  two  spheroids  in  Figs.  19a  and  19b  are  added  together 
the  result  is  the  even-parity  fluence  and  the  polar  plot  of  this  addition 
is  presented  in  Fig.  19c. 

The  selected  synthesis  function  must  have  the  general  form  repre¬ 
sented  by  Fig.  19c.  In  Chapter  II  a  centered  ellipsoid  was  proposed 
as  an  angular  synthesis  function.  Though  this  function  would  adequately 


Angular  Synthesis  Function 


approximate  the  angular  distribution  of  Fig.  19c,  it  does  not  meet  the 
first  criterion  as  already  demonstrated.  A  function  that  would  accurately 
approximate  Fig.  19c  is 

A  U(4)  =  /  +  (G.2) 


A  polar  plot  of  this  function  appears  in  Fig.  19d.  A  comparison  of  Fig. 
19c  and  Fig.  19d  demonstrates  the  adequacy  of  Eq.  (G.2)  to  approximate 
the  angular  distribution  of  the  even-parity  fluence  in  an  air-over- 

* 

ground  problem. 

The  largest  component  of  the  Boltzmann  fluence  in  an  air-over-ground 
problem  is  aligned  parallel  with  a  ray  drawn  from  the  spatial  location 
of  the  point  source  to  the  spatial  location  of  interest.  The  angular 
coordinates  of  this  streaming  direction  in  the  local  coordinate  system 
changes  from  one  spatial  location  to  another.  Let  the  streaming  direction 
for  a  spatial  location  be  represented  by /Jq  and  in  the  local  coordin¬ 
ate  system.  If  the variable  in  Eq.  (G.2)  is  defined  in  terms  of  the 


law  of  cosines,  the  following  results 


A  -  I  +  C  fJ-oU.  *  v/" i-mJ1  J  i-m1  cosC'X-'Xf,)) 


(G.3) 


Eq.  (G.3)  represents  the  angular  shape  presented  in  Fig.  19d  and 
has  the  capability  of  being  aligned  in  the  streaming  direction  at  any 
spatial  location.  Fig.  20  is  a  polar  plot  of  Eq.  (G.3). 

The  synthesis  function  represented  in  Eq.  (G.3)  satisfied  the  two 
criteria  previously  presented.  In  the  collocation  method  Eq.  (G.3)  is 
used  in  the  following  form 
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The  coefficients  and  are  determined  from  solving  the  global 
matrix  formed  from  applying  the  collocation  method  to  the  EPFBE.  .These 
two  coefficients  can  be  adjusted  to  represent  even-parity  angular  flu- 
ence  shapes  varying  from  highly  anisotropic  to  isotropic.  The  form  of 
Eq.  (G.4)  is  simple  and  was  analytically  integrable  when  placed  in  the 
EPFBE  used  in  the  collocation  method. 
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APPENDIX  H 


FLOW  CHARTS 

This  appendix  contains  the  flow  diagrams  for  both  the  Galerkin  and 
collocation  program  presented  in  Cahpters  II  and  III  respectively  for 
generating  a  matrix  representation  of  the  EPFBE.  An  explanation  of  the 
numbers  in  the  diagrams  is  presented  below. 

Galerkin  Program  (Fig.  21) 

1.  The  following  parameters  are  inputed. 

»  * 

a.  Maximum  number  of  mesh  lines  in  ^3,  2,u(,X  directions. 

b.  Number  of  spatial  and  angular  test  functions  per  finite  element. 

c.  Global  spatial  node  and  angular  node  matrix  that  relates  nodes 
to  elements. 

d.  Mesh  node  values  forp.z,^,  and  X  coordinates . 

e.  Order  of  differential  scattering  cross  section. 

f.  Cross  sections. 

2.  Quadrature  sets  for  use  with  spatial  and  angular  numerical  integration 
are  selected  and  inputed. 

3.  Print  out  all  input  data  as  a  check. 

4.  Zero  all  elements  in  the  global  matrix  and  load  vector. 

5.  Determine  coordinates  for  defining  all  angular  finite  elements. 

—  Start  loop  over  spatial  finite  elements. 

6.  Select  a  spatial  finite  element  and  define  its  location  in  the 
spatial  mesh. 

—  Start  loop  over  angular  finite  elements. 

7.  Select  an  angular  finite  element  and  from  step  5  identify  its  loca¬ 
tion  in  the  angular  mesh. 
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8.  From  the  coordinates  of  these  selected  elements  calculate  the  para¬ 
meters  needed  to  transform  each  into  a  canonical  element  representation. 

9.  Select  a  quadrature  point  for  the  ,  and  Tc  variables  and  an 

accompanying  weight  function. 

10.  Select  a  combination  of  bilinear  Lagrange  polynomials  representing 
the  spatial  and  angular  test  functions. 

11.  The  numerical  evaluation  of  the  integrals  involved  with  the  source 

term  of  the  Galerkin  weak  form  equation  is  performed.  The  integrands  of 

* 

these  integrals  are  expressed  in  terms  of  the  test  function  and  their 
derivatives  selected  in  step  10  and  are  evaluated  at  the  quadrature 
points  defined  in  step  9.  These  evaluated  terms  are  multiplied  by  the 
corresponding  weight  function  associated  with  these  quadrature  points  and 
the  product  is  stored  in  an  element  of  the  load  vector.  The  location  of 
this  element  in  the  load  vector  is  determined  by  the  test  function  com¬ 
bination.  This  element  represents  the  evaluated  integrals  in  the  source 
term  when  all  quadrature  points  and  weights  have  been  used  and  their 
resultant  products  summed.  The  value  of  the  test  functions  and  their 
derivatives  that  appear  in  the  source  term  integrands  are  determined  by 
a  series  of  function  routines  and  subroutine  programs.  These  subprograms 
define  the  test  functions  in  a  canonical  element,  evaluate  these  func¬ 
tions  for  the  selected  quadrature  points,  and  asssemble  these  evaluated 
functions  into  a  representation  of  the  source  term  integrands. 

12.  Select  a  combination  of  bilinear  Lagrange  polynomials  representing 
the  spatial  and  angular  trial  functions. 

13.  The  numerical  evaluation  of  the  local  integrals  containing  the 
even-parity  fluence  in  the  Galerkin  weak  form  equation  is  performed. 

The  integrands  of  these  integrals  are  expressed  in  terms  of  the  trial 
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functions  and  test  functions  selected  in  steps  12  and  10  respectively. 
These  integrands  are  evaluated  at  the  quadrature  points  defined  in  step 
9  and  the  result  of  this  evaluation  is  multiplied  by  the  weight  functions 
associated  with  the  quadrature  points.  This  product  which  represents  a 
partial  evaluation  of  the  weak  form  integrals  is  stored  in  one  element 
of  the  global  stiffness  matrix.  The  position  of  this  element  is  deter¬ 
mined  by  the  test  and  trial  function  combinations.  This  element  repre¬ 
sents  the  full  evaluation  of  the  Galerkin  weak  form  local  integrals  for 
the  selected  test  and  trial  function  combination  after  the  contribution 
of  all  quadrature  point  evaluations  have  been  summed.  As  mentioned  in 
step  11  several  subprograms  contribute  to  this  evaluation. 

14.  Has  a  different  spatial  trial  function  been  selected  in  step  12 
since  the  last  iteration? 

15.  Select  an  angular  finite  element  for  global  evaluation  of  the  global 
integrals  in  Galerkin  weak  form  equation. 

16.  Calculate  the  parameters  needed  to  transform  the  angular  element 
selected  in  step  15  into  a  canonical  element  representation. 

17.  Select  quadrature  points  for  the  u'  and  x'  variables  and  the  accom¬ 
panying  weight  functions. 

18.  Select  a  bilinear  Lagrange  polynomial  trial  function  defined  in  the 
angular  element  selected  in  step  15.  The  test  functions  and  spatial 
trial  functions  are  the  same  as  previously  chosen. 

19.  The  numerical  evaluation  of  the  global  integrals  (integrals  defined 
over  the  entire  angular  domain)  resulting  from  the  Ku  and  Gg  operators 
is  performed.  The  integrands  of  these  integrals  are  expressed  in  terms 
of  the  test  function  selected  in  step  10,  spatial  trial  function  selected 
in  step  12,  and  angular  trial  function  in  step  18.  The  result  of  this 


115 


evaluation  is  stored  or  added  to  existing  values  in  the  global  stiffness 
matrix  in  a  position  .determined  by  the  test  and  trial  function  combina¬ 
tion.  The  subprograms  mentioned  in  step  11  and  13  perform  the  evaluations. 

20.  Have  all  angular  trial  functions  for  the  global  integral  evaluation 
been  used? 

21.  Have  all  quadrature  points  for  the  u'  and  x'  variables  been  used 
to  evaluate  the  global  integrals? 

22.  Have  all  angular  elements  been  used  for  the  global  integral  evalua- 

* 

tions? 

23.  Have  all  spatial  and  angular  trial  function  combinations  in  step 
12  been  used? 

24.  Have  all  spatial  and  angular  test  functions  combinations  in  step 
10  been  used? 

25.  Have  all  quadrature  points  for  the  p  ,2  ,/*,  and  X  variables  been 
used? 


26.  Have  all  angular  elements  from  step  7  been  used? 

27.  Have  all  spatial  elements  from  step  6  been  used? 

Collocation  Program  (Fig,  22) 

1.  The  following  parameters  are  inputed. 

a.  Number  of  collocation  points  in  p  direction. 

b.  Number  of  collocation  points  in  Z  direction. 

c.  Number  of  mesh  points  in  p  direction. 

d.  Number  of  mesh  points  in  Z  direction. 

e.  r  coordinate  collocation  points. 

f.  Z  coordinate  collocation  points. 

g.  Mesh  coordinates  in  p  direction. 

h.  Mesh  coordinates  in  Z  direction. 
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i.  Total  macroscopic  cross  section  and  P^  cross  sections. 

j.  Height  of  burst. 

2.  All  global  matrix  elements  zeroed.  • 

3.  Start  of coordinate  collocation  loop.  A  p> coordinate  collocation 
point  is  selected  and  used  for  evaluation. 

4.  Splines  in  p  coordinate  that  have  a  value  at  selected  p  coordinate 
mesh  point  are  identified. 

5.  Start  of  2  coordinate  collocation  loop.  A  2  coordinate  collocation 
point  is  selected  and  used  for  evaluation. 

6.  The  streaming  direction  and  90  degrees  off  the  streaming  direction 
in  the  x=0  plane  at  the  selected  spatial  collocation  points  are  deter¬ 
mined. 

7.  Ellipsoidal  synthesis  function  pitch  is  set  based  on  the  streaming 
direction.  Spatial  derivatives  of  the  pitch  function  ( cos^or Ji'i'L  ) 
are  determined  for  the  selected  spatial  collocation  points. 

8.  Splines  in  2  coordinate  that  have  a  value  at  selected  2  coordinate 
are  identified. 

9.  Start  loop  over  all  splines  in  p  coordinate  that  were  previously 
identified  in  4. 

10.  Start  loop  over  all  splines  in  2  coordinate  that  were  previously 
defined  in  8. 

11.  Evaluate  all  terms  in  EPFBE  that  have  spatial  dependence.  This 
encompasses  an  evaluation  of  the  Q  ((»  Qj  c* j  terms  in  Eq.  (3.11)  and 
their  derivatives. 

12.  Decision  block  to  determine  if  all  splines  in  2  coordinate  identified 
in  8  have  been  evaluated  (Y=Yes,  N=No). 
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13.  Decision  block  to  determine  if  all  splines  in  p  coordinate  identi¬ 
fied  in  4  have  been  evaluated. 

14.  Evaluate  all  terms  in  EPFBE  that  have  angular  dependence.  This 
encompasses  an  evaluation  of  the  a.i;  +  a,,.  -  (jul/u0  * 

J  J 

term  in  Eq.  (3.11)  and  its  derivative  with  respect  to  the  u  and  x  vari¬ 
ables.  (Primed  and  unprimed.) 

15.  The  spatial  and  angular  contributions  of  the  EPFBE  found  from  11 
and  14  respectively  are  combined  to  calculate  the  value  of  all  terms  in 
the  EPFBE. 

16.  Value  of  EPFBE  calculated  in  15  is  stored  in  appropriate  global 
matrix  element.  The  value  of  the  source  terms  is  stored  in  the  load 
vector.  The  position  of  these  elements  in  the  global  matrix  and  load 
vector  is  determined  by  the  spatial  spline  combination  and  angular  col¬ 
location  points. 

17.  Decision  block  to  determine  if  both  angular  collocation  points  have 
been  evaluated. 

18.  Decision  block  to  determine  if  all  Z  coordinate  collocation  points 
have  been  evaluated. 

19.  Decision  block  to  determine  if  all  p  coordinate  collocation  points 
have  been  evaluated. 


Global  matrix  formed. 


1, 
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APPENDIX  I 

DERIVATION  OF  NORMALIZED  B. SPLINES  OF  ORDER  4 
In  Chapter  III  the  spatial  variation  of  neutral  particle  fluence 
was  approximate  by  cubic  spline  functions  defined  in  the  and  Z  variables 
The  value  of  these  functions  for  various  spatial  locations  was  calculated 
in  the  collocation  program  by  a  subroutine  called  BS.  The  following 
derivation  generates  the  equations  used  in  this  routine  for  defining 
and  evaluating  these  splines. 

Assume  the  node  sequence  t1 ,  t2,  t3>  t4,  t5  illustrated  in  Fig.  23. 


Fig.  23.  B  Spline  Node  Sequence 

<-  I 

Define  the  expression  (s- 1)+  such  that 
<  -  order  of  the  spline 
<-  <  -  degree  of  the  spline 


(S-tC 


is-i) 


if  5  2  t 

1-f 


Using  the  above  definitions  define  a  new  variable  as 

n:,kw=  . 


(i.U 


where 

i.  =  number  of  the  first  node 

CV-AJ-  divided  difference  operator  operating  on  the  S  variable  and 
defined  In  the  following  manner. 

L*Js  fc si  - 

_ a  (  C  . "f(Sj  '  CX, . (xn-xe) 


Using  Eq.  (1.1)  and  setting  K=4  the  following  representation  of  the  B 
splines  can  be  written  for  the  node  sequence  in  Fig.  23 


(1.2) 


Now  apply  the  definition  of  the  divided  difference  operator 

s  ts-t, 


(1.3) 


The  following  series  of  equations  represent  a  repeated  application  of 
the  divided  difference  operator  to  the  two  terms  on  the  right  hand  side 
of  equation  1.3 

rt.  f  t,l.  I*  +  i J_  r+  f.  4-  7  + 

(1.4) 


A  _  r.  .  j  i  i  .  '  ,  ,J  _  L  tj,ts,tr  Js  (S-t  )♦- Ctj.fj.f  t)  * 

A-  - - * - 7 - - - 

1-5  -  T* 

O  r  +  +  j.  +  rc  4\3  -  (*-*!«•  -  ts-t  )* 

'-s  -  X  , 

To  A  and  B  the  following  reductions  can  be  made 
^  r  j.  ,  ,  .  ..3  Lt*.*rJs(s-tli  +  Cti.tls  Cs-t)l 

C  =  i  - 


F  =  Lt.,u,tj]s  t s-oi  = 


ij.  t, 


(1.5) 


(1.6) 


0=  F*  •  . c -  C t. .O. u-tli —  (1 ,7) 

X*,  -  x2 


(1.8) 


Applying  the  defintion  of  the  divided  difference  operator  to  C,  D,  E, 
and  F  results  in 


G  =  Ct 

tj-  t»* 


/-i  ■  c-s-t;'  = 


tj  .  ^ 


(1.9) 

(1 .10) 
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n 


X  =  C  tj,  tjlj  (s~i  ),.  = 


J=  L±,t±Js(s-t)i-. 


(t3-t)+  - 

tj  -  tz 

(ii-t)i  -  H.-iK3 

■tiV'.t, 


(1. 11) 

(1.12) 


By  back  substituting  the  following  relationships  are  formed. 

G-  H 


C  * 


tj  -  tj 
H  -  I 


F  = 


X- J~ 

vt. 


(1.13) 

(1.14) 

(1.15) 


Further  back  substituting  results  in 

«•  -fefr- 


8  * 


D  -  F 

tu-t, 


(1.16) 


(1.17) 


Finally  the  form  of  the  B  spline  used  in  the  collocation  program  is 
defined  as 

~  A  -  8 


=  (ts-t) 


(1.18) 


The  spline  in  Eq.  (1.18)  is  normalized  such  that 

rts 

I  MIS  (tj  clt  -  I 


(1.19) 
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